4- 


END 

DATE 

riLMEO 

10-77 

DDC 


OQC  FILE  COPY  ADA  044  76 


AFGL-TR-77-0129 


SIMULATION  STUDY  OF  AIRBORNE  GRADIOMETRY 


Klaus-Peter  Schwarz 


The  Ohio  State  University 
Research  Foundation 
Columbus,  Ohio  43212 


May  1977 

Scientific  Report  No.  8 


Approved  for  public  release;  distribution  unlimited 


Am  FORCE  GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  C0MI^1AND 
UNITED  STATES  AIR  FORCE 
HANSCOM  AFB,  MASSACHUSETTS  01731 


Qualified  requestors  may  obtain  additional  copies  from  the  Defense 
Documentation  Center.  All  others  should  apply  to  the  National 
Technical  Information  Service. 


security  classification  of  this  page  ''l*?i<in  Omit  Enlerrd) 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


■?  REPORT  DOCUMENTATION  PAGE 


2.  GOVT  ACCESSION  NO 


AFGLHTR-77 


5 TYPE  OF  RtaoBT  A Ptmoa  covebeo 

Scientific  Interim  /.  ^ , 


4.  title  (snij  Subtitle) 


SIMULATION  STUDY  OF  AIRBORNE  GRADIOMETRY, 


Scientific  Report  No 


6.  PERFORMING  ORG  REPORT  NUMBER 

l"1!Dept.  Geodetic  Science  No.  I 


7.  AuTHORfJ) 


F19628 


Klaus-Peteiy  Schwarz 


PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Department  of  Geodetic  Science 
The|Ohio  State  University-1958  Neil  Avenue 


APE  4 4 UNI 

62101 F 
760(^1302 


umbus.  Ohio  43210 


CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Air  Force  Geophysics  Laboratory- 
Hanscom  AFB,  Massachusetts  01731 


14.  monitoring  agency  name  4 AODRESSfi/  trom  Controtlmi  OHicf 


Unclassified 


16.  distribution  statement  (of  this  Report) 


Approved  for  public  release;  distribution  unlimited 


7.  Distribution  statement  (of  the  abatrmet  erttereJ  in  Block  20,  if  different  from  Report) 


TECH,  OTHER 


airborne  piradiometry,  downward  continuation,  error  model  (Markon 
least-squares  collocation. 


'0.  abstract  rContInum  .in  reversm  aide  If  neceeaarv  end  iJeniify  .**v  .‘if.jrJt 


A simulated  field  of  gravity  anomalies  and  sceond-order  gradients  is 
used  to  study  the  recovery  of  gravity  anomalies  at  ground  level  from  mea- 
surements at  tlight  level. 


The  spectral  properties  of  gravimetric  quantities  simulatctl  by  an  array 
of  mass  points  arc  discussed  for  first  and  second  order  gradients  and  some 


PORM 
t JAN  73 


OITion  1 NOV  65  IS  OBSOLETE 


Unclassified 


SECURtTY  CL  ASSiFtCATlO-l  QF  THIS  F- AGgr***h»n  Datm  Enturmd) 


errors,  uncorrelated  normally  distributed  errors  and  correlated  errors 
from  a stationary  Markov  sequence  are  used.  Results  agree  well  with 
estimates  obtained  in  a previous  accuracy  study. 

Computations  show  that  a fast  operational  program  can  be  obtained 
by  using  a least- squares  collocation  procedure.  The  programs  and  a 
sample  computation  are  part  of  the  report. 


Unclassified 


security  classification  or  this  PACE^UTi..,  Dmf  Fnl»r»tl) 


1 

k 


s 


FOREWORD 

This  report  was  prepared  by  Dr.  Klaus-Peter  Schwarz,  assistant 
to  Dr,  Helmut  Moritz,  Professor,  Technical  University  at  Graz  and  Adjunct 
Professor,  Department  of  Geodetic  Science  of  The  Ohio  State  University,  under 
Air  Force  Contract  No.  F19628-76-C-0010,  The  Ohio  State  University  Research 
Foundation  Project  No.  4214A1,  Project  Supervisor,  Urho  A.  Uotila,  Professor, 
Department  of  Geodetic  Science,  The  contract  covering  this  research  is  ad- 
ministered by  the  Air  Force  Geophysics  Laboratory  (AFGL),  Hanscom  Air 
Force  Base,  Massachusetts,  with  Mr.  Bela  Szabo,  Project  Scientist. 


I 


CONTENTS 

1.  Introduction 

2.  Mass  Models  and  Their  Spectral  Properties 

3 . Error  Model s 

4 . Results 

5.  Brief  Description  of  Programs 

6.  Conclusions 
References 

Appendix  A : Computer  Programs 
Appendix  B : Sample  Computations 


i'.' 


1.  Introduction 


In  a previous  report  (Schwarz,  1976  ) the  accur-acy  of 
airborne  gradionetry  has  been  studied  and  some  conclusions  have 
been  drawn  about  optimal  point  configurations  and  data  combi- 
nations. This  report  supplements  some  of  the  previous  investi- 
gations. Not  much  can  be  added  with  respect  to  the  expected 
accuracy.  Simulation  studies  display  the  behaviour  of  individual 
experiments  only  and  are  therefore  not  suited  to  check  results 
of  an  accuracy  study.  The  interest  of  a simulation  study  is  there 
fore  not  so  much  in  the  field  of  accuracy  but  in  the  domain  of 
operational  realization  and  optimal  performance. 

In  order  to  get  an  operational  program  for  airborne 
gradiometry  the  most  important  problem  to  cope  with  is  the 
efficient  handling  of  large  amounts  of  data.  The  proposed 
measuring  system  will  produce  about  250  observations  per  profile 
and  degree.  In  order  to  cover  a 20°  x 25°  area  with  profiles 
spaced  at  1°  we  have  to  treat  130  000  measurements.  For  mean 
gravity  values  below  1°  x 1°  we  have  to  use  20'  spacings  and 
the  above  number  of  measurements  will  triple.  It  has  been  shown 
in  Schwarz  (1976)  how  the  number  of  observations  can  be  reduced 
without  significantly  imoairing  the  accuracy  of  the  results.  For 
an  operational  program,  however,  it  will  be  necessary  to  use  all 
information  available.  Not  so  much  to  increase  accuracy  but  to 
make  results  more  reliable.  Therefore,  we  have  to  incorporate 
both  viewpoints  in  such  a program.  Efficiency  asks  for  data 
selection  in  each  computation  step.  Reliability  requires  the 
processing  of  all  data  available. 

A second  reason  to  carry  out  a simulation  study  is  the 
treatment  of  pathological  error  situations.  Accuracy  studies  are 
usually  oerformed  under  the  assumption  that  the  cbs erv a t i ona 1 
errors  follow  a Gaussian  distribution.  '/Jith  complex  measuring 
systems  on  a moving  base  this  assumption  may  not  be  realistic. 
Besides  correlated  errors  biases  a''e  likely  to  occur  in  air  some 
gradiometry.  Effects  of  this  kind  are  easy  to  simulate  and 
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results  are  important  for  the  planning  of  experiments.  Thus, 
the  actual  development  of  the  error  budget  coming  from  different 
sources  will  help  to  plan  an  effective  updating  procedure. 

Finally,  the  checking  of  such  a program  in  a controlled 
experiment  is  a worthwile  exercise  by  itself.  Not  only  because 
of  its  complexity  but  because  instabilities  stemming  from  the 
downward  continuation  problem  should  be  controlled  in  the  best 
way  possible.  Special  care  must  be  taken  that  the  methods  for 
generating  the  data  are  truly  independent  of  the  data  processing 
procedures.  If  this  can  be  achieved  simulation  studies  will  give 
a reliable  base  to  handle  real  data.  Large  differences  in  the 
results  from  actual  and  from  simulated  data  will  indicate  that 
the  mathematical  model  needs  refinement.  The  nature  of  the  re- 
finement can  often  been  guessed  from  the  simulation. 

The  numerical  treatment  of  the  problem  requires  the 
consideration  of  the  following  four  steps: 

Step  1 ...  Gravity  and  gradiometry  data  are  generated  at  ground 
and  at  flight  level. 

Step  2 ...  Data  at  flight  level  are  corrupted  by  the  error 
model  . 

Step  3 ...  Gravity  anomalies  at  ground  level  are  estimated  from 
data  at  flight  1 evel . 

Step  4 ...  Estimated  and  model  anomalies  are  compared  at  ground 
level . 

These  steps  will  be  covered  in  the  next  three  sections. 
Data  simulation  in  section  2,  error  models  in  section  3,  and 
estimation  and  comparison  in  section  4.  Section  5 will  give  a 
short  review  of  the  programs  which  are  listed  in  an  appendix. 

2.  Mass  Models  and  Their  Spectral  Properties. 

The  basic  assumption  underlying  the  simulation  of 
gravimetric  quantities  by  a point  mass  model  can  be  formulated 
in  the  following  way:  The  field  generated  by  such  a model  in  a 
limited  region  can  be  regarded  as  a sufficient  approximation  o^ 
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the  anomalous  gravity  field  in  this  region.  In  a number  of 
applications  the  actual  field  can  only  be  described  by  statis- 
tical parameters.  In  such  a case  the  above  would  imply  that  the 
statistical  properties  of  the  simulated  field  can  also  be  con- 
sidered as  approximations  of  the  actual  quantities.  This  as- 
sumption has  important  theoretical  and  practical  implications. 
Some  of  them  will  be  examined  in  the  sequel. 

Let  us  first  consider  the  simple  case  that  all  anomalous 
masses  are  concentrated  on  a plane  at  depth  and  that  the 

simulated  function  is  wanted  on  a parallel  plane  . Using 

Cartesian  coordinates  we  can  write 


— 00 

where  density  function  representing  the  anomalous 

masses  s the  upward  continuation  operator,  and  g is  the 

simulated  gravity  function.  Since  we  will  use  the  spectral 
density  function  later  on  we  will  call  f the  mass  density 
function  in  the  sequel.  Obviously,  g will  be  singular  at  z^  , 
i.e.  it  is  only  defined  above  this  plane,  and  we  have  the 
conditon  z,  > z^^  . 

Equation  (2.1)  defines  a double  convolution  of  the 
functions  h and  f which  may  be  written  as 

g(x,y,z_.)  = h(x,y,Z2-z^)-»f(x,y,zp  (2.2) 

where  » is  the  convolution  symbol.  Since  we  consider  only  inte- 
grations in  the  (x,y)-plane  we  will  write 

g(x,y)  = h(x  ,y)--f (x  ,y)  , (2.3) 

keeping  in  mind  that  both  g and  f refer  to  a fixed  z . Let 
us  assume  that  the  Fourier  transforms  of  all  three  functions  in 
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equation  (2.3)  exist.  We  will  denote  them  by  capital  letters 
and  reserve  small  letters  for  the  data  domain.  As  an  example  we 
have 


G(u,v)  = /.'g(x,y)e  ^ dxdy 


as  Fourier  transform  of  g(x,y)  and  the  inverse  relation 


9(x,y)  = — //G  ( u , V ) e^  dudv  . 

4 T - OO 


G(u,v)  is  also  called  the  spectrum  of  g(x,y)  . Forming  all 
three  transforms  we  can  make  use  of  the  simple  relations  in  the 
spectral  domain 

G(u,v)=H(u,v).F(u,v)  (2.6; 

and  using  the  inversion  formula  (2.5)  we  obtain 


g(x,y)  = -^  //H(u,v)F(u,v)e^^'"'''"''''’’dudv  . 
4 TT 


The  transfer  function  H(u,v)  is  determined  by  the  geometrical 
relations  between  the  planes  and  and  is  of  the  form 


+ aj.  -i'u(x-.'c')+v(y-y')} 

H(u,v)  = l\— ^ ^ VS,-  dx'dy 

-oo{  (x-x-)‘-+(y-y')‘-  + d‘-'^^^“ 


where  d=z,  -z^. 

Evaluating  the  integral  we  obtain 


H(u  ,v)  = 2-e 


-av  u + V 


(2.9) 


and  we  can  write  formula  (2-7)  as 


9(x,y)  = 


I \/~2  F 

1 . -.-d/u  -V  p^^^^^gi(ux^vy) 


dud  V 


(2.10) 


Equation  (2.10)  shows  that  the  transfer  function  is  a smoothing 
function  which  affects  the  high  frequencies  most.  The  degree  of 
the  smoothing  is  dependent  on  the  size  of  d , i.e.  on  the 

separation  of  the  two  planes.  To  illustrate  this  point  table  2.1 

gives  smoothing  factors  for  different  values  of  d and  different 
frquencies.  To  make  the  results  applicable  to  the  mass  model 
used  later  on,  a grid  of  61  by  61  equidistant  mass  points  has 
been  chosen.  The  variable  d is  expressed  in  units  of  the  mass 

point  spacing.  Since  H(u,v)  has  circular  symmetry  we  have  used 


H ( w ) = 2 1 e 


(2.11) 


where 


w =\ 

w 

d 

0.5  1 

1.0 

1 . 5 

2.0 

2 . 5 

1 

.950  j 

.902 

.357 

.814 

.773 

5 

.773 

. 597 

. 462 

.357 

. 276 

10 

. 597 

. 357 

.213 

.127 

.076 

15 

. 462  ' 

i 

.213 

.099 

.045 

. 021 

20 

1 

.3  57  1 

.127 

. 045 

.016 

.005 

25 

.276 

. 076 

. 021 

.006 

.002 

30 

.213 

. 045 

.010 

.002 

. 000 

Table  2.1  Smoothing  of  gravity  anomaly  soectrum 


It  is  apparent  from  this  table  that  most  of  tne  high 
frequency  information  is  lost  if  d becomes  larger  than  the 


mass  point  spacing.  This  will  result  in  a very  smooth  gravity 
field  because  its  structure  is  determined  by  a few  low  fre- 
quencies only, 

A similar  consideration  applies  for  second-order 
derivatives.  We  will  show  it  for  the  vertical  derivative  of 
Using  equation  ^2,10)  we  obtain 


^ ^ 1 ^ ' 


F(u,v)e"‘^^^'^^’dud 


Thus 


H,  (u  , v)  = l/u^  + v^  H(u  , v' 


H (w)  = w H(w) 


Results  are  given  in  table  2.2. 


ini 


0.5 

1 . 0 

h 

1.5 

2.0 

2.5 

. 950 

. 902 

1 1 

.857  ! 

.314  i 

1 

.773 

3.865 

1 

2.987 

2.309 

1.785  i 

1.380 

1 

5.975 

3.570 

2.133 

I 1.274 

j 

.761 

6.923 

1 

3.200 

1.478 

.682 

.315 

7.140  i 

1 

2.549 

1 

.910 

. 325  ! 

1 

.116 

6.399  i 

1.904 

.525 

! .145  1 

.040 

6.399 

1.365 

1 

.291 

: .062  1 

.013 

Tabel  2.2  Smoothing  of  second-order  vertical 
gradient  spectrum 


7 


Since  the  high  frequencies  of  the  second-order  gradients 
are  strongly  amplified  the  smoothing  effect  of  H(w)  is  counter- 
balanced if  d does  not  become  too  large.  Thus,  even  with  d 
twice  the  mass  point  spacing  we  can  expect  adequate  information 
on  the  frequencies  up  to  15  or  20.  There  is  not  enough  empirical 
information  at  the  moment  on  the  actual  variation  of  the  second- 
order  gradients  to  decide  whether  a field  with  such  frequency 
content  is  too  smooth.  It  is  obvious,  however,  from  table  2.2 
that  by  varying  the  grid  density  of  the  mass  points  and  their 
depth  we  can  model  a wide  range  of  different  fields. 

So  far  no  assumptions  have  been  made  on  the  mass  density 
function  f(x',yjz^)  except  that  it  should  possess  a Fourier 
transform.  Thus,  we  can  model  e.g.  a given  gravity  anomaly  field 
in  a completely  deterministic  way  by  properly  distributing  point 
masses  at  certain  depths.  However , if  the  field  to  be  simulated 
can  only  be  characterized  by  statistical  parameters  a different 
approach  must  be  taken.  In  such  a case  we  are  looking  for  a mass 
distribution  which  will  generate  a field  with  the  desired  statisti- 
cal properties.  Obviously,  the  stochastic  characteristics  of  the 
mass  density  function  and  the  influence  of  the  transfer  function 
must  be  taken  into  account. 

Let  us  consider  a two-dimensional  wide-sense  stationary 
process  f(x,y)  , i.e.  a process  which  has  constant  mean  value 
and  an  autocorrelation  r which  depends  only  on  •:  = x^-x, 
and  - = y. -y.,  . Thus,  it  is  character  i zed  by  its  first  and  second 

i 2 

moments 


E{f(x,y)}  = const. 

E'f(x  + ;,y  + -)  • f(x,y)}  = (2.  IT) 


where  E is  the  statistical  expectation  and  the  over bar  denotes 
the  complex  conjugate.  In  the  sequel  we  will  only  use  processes 
with  mean  values  equal  to  zero.  Therefore,  no  distinction  is 
necessary  between  the  autocorrelation  r_._.  and  the  autocovariance 
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Cf.  and  we  will  always  use  the  latter  one.  The  Fourier  transform 
of  is  given  by 


Sff 


(u,v) 


= IJCff  (C  ,n  )e' 


■ i ( uC  +vn ) 


dc  dn 


(2.15) 


where  S^.  is  called  the  spectral  density.  It  can  be  shown  that 
the  spectral  density  of  an  arbitrary  process  is  nonnegative 


S (u  ,v)  > 0 . 

r 1. 


In  general,  a stationary  process  f(x,y)  does  not  have  a spectral 
representati on  of  the  form  (2.5).  Thus,  equation  (2.15)  must  be 
regarded  as  the  basic  spectral  relation  for  a stationary  process. 
Given  a positive  function  S^^(u,v)  or,  equivalently,  a positive 
definite  function  c.^(^,n)  . we  can  find  a stochastic  process 
having  S^^(u,v)  as  spectral  density  and  c.^(C,n)  as  co- 
variance  function.  If  the  covari ance’ functi on  does  not  have  a 
Fourier  transform  it  can  usually  be  represented  by  a Fourier- 
Stieltjes  integral  and  the  inversion  of  thi  s ■ i ntegra  1 is  possible 
by  generalized  transform  methods.  Incidently,  this  is  the  method 
also  used  for  the  spectral  representat i on  of  the  process  itself. 

Let  us  now  consider  the  convolution  (2.3).  It  would  be 
quite  advantageous  to  have  g(x,y)  stationary  because  this  would 
allow  us  to  characterize  the  simulated  field  by  its  first  two 
moments.  One  condition  which  secures  stationarity  of  g(x,y)  is  the 
following;  If  the  process  f(x,y)  is  wide-sense  stationary  then  the 
output  of  the  convolution  (2.2)  will  also  be  stationary.  Further- 
more, f(x,y)  and  g(x,y)  will  be  jointly  stationary,  i.e.  the 
joint  statistics  of  f(x,y)  and  g(x,y)  will  be  the  same  as  the 
joint  statistics  of  f(x  + ;,y  + '-. ) and  g(x  + -;,y  + ")  . Thus  we  have 


E{  f (X  + -;  ,y  + n)g(x,y) } 


9 


or 


Similarly 


= c (;,n)»*h(5,n) 
gg  r g 


(2.16) 


(2.17) 


Forming  the  double  transforms 

Cfj('^>")>  *^,Tg(^.!n)j  c^^(:i,n)  . 

(2^7) 


S-^(li,v),  S_.^(u,v),  S (u,v)  of 

r - j.  g 

we  obtain  from  (2.6),  (2.16)  and 


S 

gg 


S 

gg 


(U.V) 
(U.V) 
(u  , v) 


S_-.-(u,v)  H(U  ,v) 

S-_(u,y)  H(u,v) 
i^g 

S.^-(u,v)  H(u,v),^  . 

X.  4. 


By  use  of  equation  (2.9) 


S ( u , 7 
? ■? 


= S,,(u,v) 


(2.18) 


(2.19) 


Formula  (2.19)  shows  that  the  spectral  density  of  the  simulated 
process  is  related  in  a simole  way  to  the  spectral  density  of  the 
mass  density  function  if  f(x,y)  is  stationary.  The  importance 
of  this  result  has  already  been  stressed  by  Naidu  (1968)  and 
Grafarend  (1970).  It  is  only  valid  for  the  planar  approximation 
while  for  the  spherical  case  a correction  term  is  necessary. 

Using  the  inverse  we  obtain 


c 

gg 


(-,n) 


( u , V ) e" 


' dudv 


(2.20) 


and  for  the  variance 
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c (0,0)  = //e-^^vu  s (u,v)dudv  . (2.21) 

y y _oo 

There  are  simple  relations  between  the  spectral  densities  of 
different  first  and  second  order  gradients  which  are  e.g.  derived 
in  Naidu  (1968)  and  Kubackova  (1974). 

In  order  to  select  a special  model  we  have  now  to  choose 
a specific  spectral  density  or,  equivalently,  a specific  covari- 
ance function.  A number  of  different  Gaussian  models  have  been 
considered  by  Grafarend  (1970).  In  case  of  point  mass  anomalies 
a model  correspond i ng  to  white  noise  seems  to  be  most  adequate. 

In  two  dimensions  it  is  sometimes  called  an  incoherent  process 
and  is  defined  by  the  covariance  function 

Cff (Xi  ,yi  ,y2)  = q ( x ^ ,y ^ ) 5 ( x, -x ^ ) 5 (y^ -y ^ ) (2.22) 

where  5 denotes  the  0 i rac-f unct f on , 

To  get  a stationary  process  we  must  have 

q(x,y)  > 0 . (2.23) 

Since  white  noise  processes  have  infinite  intensity 


E{if(x,y).^}  = - 

it  is  useful  to  define  the  average  intensity  of  the  output  which 
in  our  case  is 

E(  ' 9(x,y)  ; “I  = q(x,y)»»|  h(x,y)  , (2.24) 

and  if  h(x,y)  takes  significant  values  only  in  a limited  region 
A near  the  origin,  we  get 

E{  ' g(x,y)  = q(x,y)j'/| 

A 


h ( X ,y ) “dxdy 


(2  . 25i 
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Another  model  which  might  be  worthwhile  considering  in  our  case 
is  a stationary  Markov  sequence. 

So  far,  we  have  used  only  one  plane  of  generating  masses. 
As  is  apparent  from  tables  2.1  and  2.2  the  combination  of  different 
planes  would  be  advantageous  for  a realistic  model  of  first  and 
second  order  gradients.  Schwahn  (1975)  has  investigated  this  case 
and  has  shown  that  the  autocovariance  function  of  the  simulated 
field  can  be  determined  by  adding  the  autocovariance  functions 
of  the  partial  fields  if  no  crosscovari ances  between  the  planes 
are  present.  This  case  seems  to  be  applicable  in  all  kinds  of 
model  computations.  If  crosscovariances  in  z-direction  cannot 
be  neglected  the  full  covariance  matrix  of  the  partial  fields  must 
be  used  and  computations  become  extremely  laborious.  The  case  of 
a thick  sheet  of  random  masses  and,  as  a special  case,  of  a semi- 
infinite medium  has  been  treated  by  iNaidu  ( 1968). 

The  practical  procedures  of  generating  gravimetric  data 
from  a point  mass  model  can  only  approximate  the  discussed  models. 
Usually  the  anomalous  masses  are  allocated  to  the  intersections  of 
a grid  in  the  (x,y)-plane,  i.e.  the  generating  field  is  a two- 
dimensional  array  of  mass  points  equidistant  along  the  axes.  If 
the  grid  covers  the  infinite  plane  a good  approximation  of  the 
above  models  is  always  possible.  If  the  grid  is  only  given  in  a 
finite  region,  the  situation  becomes  more  difficult  because  theo- 
retically we  will  loose  stationarity  of  the  simulated  function. 

It  can  be  expected,  however,  that  stationarity  will  be  good 
enough  for  all  practical  purposes  if  the  array  of  mass  points  is 
chosen  prooerly.  Obviously,  the  density  of  the  grid,  its  extension, 
and  the  depth  of  the  generating  masses  will  be  important  oara Pi- 
eters. 

3.  Error  ‘-lodels. 


Two  types  of  statistical  error  models  have  been  used  in 
the  ccmpu  t a t i on  s : normal  models  and  '-'arkov  models.  The  term 
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normal  model  will  refer  to  a series  of  numbers  taken  from  a normal 

~>  2 
distribution  Mfuja"}  with  mean  ^ and  variance  j . The  term 

Markov  model  will  be  used  in  connection  with  Markov  sequences  of 
first  and  second  order.  The  following  discussion  will  be  main- 
ly directed  towards  Markov  sequences. 

We  will  assume  wide-sense  stationarity  for  these  models 
and  thus  be  able  to  use  some  of  the  results  of  section  2.  Matters 
are  simplified  to  a certain  extent  by  the  fact  that  we  are  con- 
sidering sequences  only.  They  can  be  viewed  as  stochastic  processes 
d(t)  where  t can  take  integral  values  only.  Furthermore,  since 
we  do  not  assume  error  correlations  between  profiles  we  will  only 
consider  the  one -d i men s i ona 1 case. 

The  difference  between  the  normal  and  the  Markov  model 
lies  in  their  correlation  characteristics.  Normal  deviates  are  un- 
correlated, elements  of  Markov  sequences  are  not.  If  e^  is  an 
element  of  the  first  model,  we  have 


E{(e.-u^)(e.-u^)l  = a- 


1 j , 


(3.1) 


where  E is  again  the  statistical  expectation. 

If  d^  = d(t)  is  an  element  of  the  second  model,  we  have 

E-  -•,,)'  = (3.2) 

^ = "x  k = 1,  2,  3,  ...  , 


where  c.  is  the  k-th  autocovariance  with  the  corresponding 
correlation 


(3.3) 


J 
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Stationarity  of  the  models  is  apparent  from  the  first  two  equa- 
tions of  (3.1)  and  (3.2).  Their  difference  is  obvious  from  the 
last  equation  in  each  group.  While  any  element  e^  is  inde- 
pendent of  any  other  element  e^  , each  element  d^  of  the 
Markov  model  does  depend  on  one  or  more  of  the  previous  elements; 

all  d.  have  the  same  univariate  distribution  but  are  correlated 
1 

The  following  presentation  of  Markov  sequences  will  only 
cover  those  characteristics  which  are  interesting  for  the  sub- 
sequent computations.  Since  the  presentation  will  be  rather 
heuristic  it  may  be  necessary  to  consult  a more  detailed  expo- 
sition. The  introductory  texts  of  Parzen  (1962),  Yaglom  (1973), 
and  Kendall  (1976)  have  been  found  especially  useful. 

Markov  sequences  are  distinguished  by  their  order.  In 
an  intuitive  way  the  order  of  the  sequence  is  equivalent  to  the 
smallest  number  of  independent  values  necessary  to  describe  the 
correlation  in  the  series.  Thus,  a Markov  sequence  of  first 
order  is  of  the  form 


d . = rd . , + e . 
1 1-1  1 


(3.4) 


where  we  have  only  one  coefficient  r describing  the  correlation 
To  obtain  the  correlation  between  d^  and  d,_,  we  apoly  the 


above  formula  to  d^  ^ and  obtain 


d , . = r d 


Inserting  this  into  equation  (3.4)  we  get 


d =r"d  +re  +e 
1 i-_  1-1  1 


Continuing  in  this  way  we  see  that  the  correlation  between  o 
•and  d , is  r'"‘  and  that  we  can  write 

1 -t- 


(3.6;. 
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Keeping  in  mind  that  r is  real  ,r,  < 1 and  that  the  correlation 
function  is  symmetric  we  obtain  the  covariances  Cj^  by 


c 


k 


C 


k 


1=0 


r r 


c +k 


1-r" 


(3.7) 


2 

Thus  we  can  express  the  variance  a of  the  Markov  sequence  of 

‘‘  -> 

first  order  in  terms  of  the  variance  by 

N 


2 


(3.8) 


2 

This  shows  that  if  'r,  is  close  to  one  j,.  will  be  much  larger 

, ' M 

than  j"  . Thus  a high  correlation  of  sucessive  values  will  pro- 
duce large  amplitudes  even  if  the  disturbances  are  small. 

The  relation  between  the  spectral  density  and  the  co- 
variance  function  of  a stationary  Markov  sequence  can  be  re- 
presented in  a form  similar  to  equation  (2.15) 


= y 


• 1 w< 


(3.9) 


where  S^^(w)  is  the  spectral  density,  c^^(k)  the  covariance 
function, and  where  k can  take  integer  values  only. 
Conversely,  we  obtain 


c 


d'i 


(k) 


S,,(w) 


e'-''""dw 


(3.10) 


Using  equation  (3.7)  we  have  for  the  Markov  sequence  of  first  order 
c_j,(  t)  = c,r 


(3.11) 


where 


c 

o 


Since  the  covariance  function  is  symmetric  in  k 


we  obtain 


"jd'"*  ' '=0'" 


(3.1 


and  using  formula  (3.9) 


dd 


(w) 


i w k 


CO 


(3.13 


After  some  rearrangement  the  summation  of  the  series  will  result 
i n 


(3.14 


or 


S 


J.d 


» 


: 1 

^(l-2rcosw+r') 


(3.15 


Thus,  we  have  determined  the  spectral  density  from  the  covarianc 
f u n c t i 0 n . 


ro 
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A Markov  sequence  of  second  order,  also  called  com- 
pound Markov  sequence  or  Yule  sequence,  is  characterized  by  two 
correlations  r^  and  r ,,  . We  can  define  it  by  the  following 
equation 


d.  =-a,d.  , -a,d.  ^+e, 

1 li-l 


(3.16) 


where  the  coefficients  a^  and  a,  are  real  numbers  and  satisfy 
a^;  < 1,  ,a^!  < 1 . They  are  connected  to  the  correlations  r, 
and  r,^  by 


r^(l-r J 


1 - r 


(3.17) 


a,  = - 


1 - n: 


or  conversely 


1 + a 


a ,,  + 


(3.18) 


1+a 


The  spectral  density  of  such  a sequence  is  given  by 


I W 


(3.19) 


e - a j e -a 


Using 


t w 


-a^  - = (e^’"'-a.  )(e’^'"'-a  J 
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and  replacing  e"  by  z we  obtai 


^ ( w ) = c . 


(z-a,  ) { 1-zaj^)  (z-a, } ( 1-za, ) 


and  aoplying  partial  fractions 


S-,.(w)  = y 


( a,  -a , ) ( 1-a, a J 


I 1-a 


^ ^ + -r-i — i 


z - a , 


• a . z 


1-at  ^"-^2  ^ 


(3.20) 


The  functions  l/(l-a^z)  and  l/(l-a.,z)  are  regular  in  the  unit 
circle.  Hence  the  series  expansions  contain  only  nonnegative 
powers  of  z . The  functions  a^/(z-a^  and  a^(z-a^)  are  regul 
outside  the  unit  circle  and  have  series  expansions 


a 1' 


; -a 


n , - n 


n = 1 


for  z ^ • E^xpanding  in  this  way  and  using  equation  (3.10) 
we  can  determine  the  covariances 


c,,(k)  = 


, a , - a , ) ( 1 - a . a , ) 


- a, 


a . 


1 - a 7 


1 - a ' 


(3.21) 


The  variance  of  the  second-order  "arkov  sequence 
can  be  expressed  in  terms  of  the  variance  -,T  . We  use  equation 
(3.16)  in  tne  form 


e.  =d.  +a,d_  , +a,d_  , 
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and  take  statistical  expectations  on  both  sides 


E{e^,e.}  - ,+a^d^_^+a,d  . _ J , ( d ^ + a ^ d . _ ^ +a  _,d  . _ J } 


Substituting  and  in  terms  of  a^  and  a^  results  in 


1 + 3 . 


J = J • 
MM  N 


(l-a^){  (l  + a^)“-ap 


(3.22) 


Equation  (3.22)  can  be  used  to  determine  the  factor  c in 
f ormu la  (3.21). 

Figures  of  covariance  functions  of  Markov  sequences  and 
of  their  corresponding  spectral  densities  can  e.g.  be  found  in 
Kendall  (1976). 

Markov  sequences  of  higher  order  can  be  generated  in  an 
analogous  way.  The  basic  equation 

d.  = - a.d.  - a.d,  , -da.  + e,  (3.23) 

1 1 1-1  - 1--  n i-n  1 ' 

has  the  spectral  density 


= c 


1 w , 

e -a 


e""'-a 


(3.24) 


which  can  again  be  used  to  determine  the  covariance  function. 

The  use  of  stationary  Markov  sequences  is  greatly 
facilitated  by  the  simple  way  in  which  they  can  be  generated. 
Equations  (3.4)  and  (3.16)  can  directly  be  used  for  this  purpose. 
Subroutines  for  normal  deviates  e^  are  usually  available  in 
program  libraries.  Otherwise,  effective  methods  to  compute  normal 
deviates  are  e.g.  given  in  Hamming  (1962).  We  then  assume  d 

O 
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and  previous  terms  to  be  zero  and  run  the  series  for  a number 
of  terms  until  the  effect  of  these  initial  assumptions  has  be- 
come negligible.  From  this  point  onwards  the  series  can  be  re- 
garded as  a Markov  sequence.  To  reach  this  point  in  a small 
number  of  steps  the  normal  deviates  e.  should  be  multiplied 
" O'"  ( ^ ^ ^ to  obtain  the  right  size  of  the 

variance  for  the  values  d.  . These  ratios  can  be  obtained  from 
formulas  (3.3)  and  (3.22). 

4 . Results 

The  data  used  in  the  following  computations  have  been 
generated  by  DMA  (Howard,  1976).  It  had  been  requested  that  the 
variance  of  the  gravity  anomalies  at  ground  level  should  be 
close  to  C = 1650  mgal'^  and  that  the  variance  of  the  horizontal 
derivatives  of  Ag  at  the  same  level  should  be  about  = 100  E“. 
The  value  of  C presupposes  that  a regional  part  of  the  aravity 
field  corresponding  to  an  expansion  of  about  degree  and  order 
10  has  been  subtracted.  The  field  has  been  generated  using  a 
grid  of  mass  points  with  a spacing  of  10'  at  a depth  of  40  km.  The 
total  area  is  10°xl5'^  , i.e.  about  4150  mass  points  have  been 
used.  The  mass  points  ha'.'e  either  a positive  or  a negative  mass 
of  constant  size  or  a zero  mass.  The  selection  was  made  acco!'- 
ding  to  a normal  distribution.  More  details  on  the  simulation  cf 
the  data  can  be  found  in  Howard  (1976). 

In  view  of  section  2 this  model  has  some  res t r i c t i on s . 

The  deoth  of  the  generating  masses  is  about  2.2  times  the  mass 
point  spacing.  Judging  from  tables  2.1  and  2.2  the  high  fre- 
quency part  of  the  first  and  second  order  gradients  will  be  very 
small  and  the  respective  fields  will  be  smooth.  Furthermore . the 
a porox i ma t i on  of  stationary  white  noise  by  a jump  process  with 
three  possible  states  (+,-.0)  seems  somewhat  inadequate.  Future 
models  should  at  least  have  two  generating  planes,  one  at  a deptn 
of  about  40-60  km,  the  other  at  about  15-25  km.  The  dense  point 
spacing  is  only  necessary  on  the  uoper  plane,  the  lower  one  may 
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have  a much  wider  spacing.  The  point  masses  should  be  taken  from 
a normal  distribution  with  the  larger  variance  on  the  lower  plane. 

The  following  results  should  be  seen  with  these  re- 
servations in  mind.  Thus,  statements  about  accuracy  may  need 
some  qualification  while  the  conclusions  about  operational  perfor- 
mance should  be  fairly  general.  Since  the  latter  was  the  main 
objective  of  this  study  there  was  no  immediate  need  for  a more 
el aborate  model . 

Fig.  4.1  shows  the  simulated  anomalous  gravity  field  at 
ground  level  with  20-mgal  contour  lines.  The  field  shows  large 
variations  with  the  e.xtreme  points  at  about  plus  and  minus  120  mgal 
and  U = 1741  mgal^  . The  variance  H is  at  about  96  E"  . 

o o 

The  formula  used  to  estimate  gravity  anomalies  at  ground 
level  from  the  simulated  data  at  flight  level  is 


s 


C 

sx 


C 


- 1 


XX 


X 


(4.1) 


where  s is  the  vector  of  gravity  anomalies,  called  the  signal, 
and  X is  the  vector  of  simulated  data,  called  observations  in 
the  sequel.  It  should  be  noted  that  contrary  to  section  2 and  3 
covariance  matrices  are  denoted  by  capital  letters  to  simplify 
comoarison  with  other  publications.  C is  the  autocovariance 
matrix  of  the  observations.  Since 


t = X - n 


(4.2) 


we  have 

C = C + C 

XX  1 1 n n 


(a. 3) 


where  t refers  to  any  first  or  second  order  grad  i ent  used  as 
observation  and  where  n is  the  error  given  by  one  of  the  error 
models.  C contains  the  crosscovariances  between  signal  and 

S X 

observation,  and  since  s and  n are  considered  to  be  uncor- 
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I 
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related  we  have 
C = C 

SX  St 

Both,  s and  t are  quantities  of  the  anomalous  field.  Their 
mathematical  relation  can  be  expressed  by  integral  equations  or 
infinite  series.  In  formula  (4,1)  these  relations  are  contained 
in  the  covariance  matrices,  i.e.  all  infinite  operations  have 
been  performed  on  the  covariances  without  approximation.  Thus, 
we  have  a consistent  model  for  heterogeneous  observations  and 
this  characteristic  property  is  preserved  even  if  the  covariance 
function  is  not  optimal. 

The  covariance  function  used  in  this  report  has  been 
described  in  Schwarz  (1976).  It  has  the  advantage  of  numerical 
simplicity  and  it  agrees  well  with  statistical  estimates  of  the 
anomalous  field.  Three  of  these  estimates  have  been  used  as 
essential  parameters  for  the  covariance  function,  namely 

= 1500  mgal^  ; = 61  km  = 111  , (4.4) 

where  C is  the  variance  of  the  gravity  anomalies,  ; is  the 

O 

correlation  length  of  the  corresponding  covariance  curve,  and 
G is  the  variance  of  the  horizontal  derivatives  of  the  gravity 

O 

disturbance.  All  quantities  refer  to  ground  level.  Fig.  4.2  shows 
this  covariance  function  as  heavy  line.  The  dashed  line  refers 
to  a covariance  function  directly  derived  from  the  simulated 
gravity  field  using  a set  of  9600  points  spaced  at  715  intervals. 
Its  essential  parameters  are 

Z = 1741  mgal-  T = 52  km  Z = 96  E“  . (4.5) 

O 3 


It  will  be  used  later  on  to  determine  the  influence  of  wrong 
assumptions  in  the  covariance  function  on  the  estimation. 


Fig.  4.2  Gravity  anomaly  covariance  functions  used  in  the 
computati ons . 

The  observations  at  flight  level  (10  km)  have  been  gener- 
ated for  a flight  speed  of  500  knots  and  an  integration  interval 
of  10  sec.  Thus,  the  separation  of  data  points  is  a constant 
2.6  kn . The  profiles  are  parallel  and  run  in  east-west  direction. 
They  cover  an  area  of  10^x15^.  Generally,  the  profile  spacing  is 
1*^ , except  in  a strip  of  2‘^xl5"  where  it  is  20'.  The  choice  of 
the  east-west  profiles  was  suggested  by  the  favourable  error  be- 
haviour in  this  direction  (Meissl,  1970)  and  because  simulations 
are  especially  simple  with  such  an  arrangement.  There  are  no  basic 
changes,  however,  if  an  arbitrary  direction  is  used  as  long  as  the 
profiles  remain  parallel  and  the  data  points  have  a constant 
separation. 

It  has  been  shown  in  Schwarc  (1976)  that  the  accuracy  of 
the  estimation  is  mainly  dependent  on  the  point  configuration 
and  net  so  much  on  the  number  of  observations  used.  On  the  other- 
hand,  it  has  been  pointed  out  tnat  for  reasons  of  reliability 
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all  available  data  should  be  used.  To  combine  both  viewpoints  the 
following  approach  has  been  taken.  An  optimal  point  configuration 
is  chosen  depending  on  the  profile  spacing  and  on  the  size  of 
the  mean  anomalies  to  be  estimated.  With  the  above  conditions 
on  parallelism  of  the  profiles  and  constant  point  separation  we 
can  move  along  the  profiles  without  changing  the  operator  R 


R 


C 

SX  XX 


(4.6) 


Moving  step  by  step  from  one  set  of  data  points  to  the  next  we 
use  all  the  infromation  available.  In  each  step  we  estimate  one 
or  several  gravity  anomalies  at  ground  level  simply  by  multi- 
plying a set  of  observations  by  the  predetermi ned  matrix  R . 
When  advancing  along  the  profiles  we  get  a whole  series  of  esti- 
mation points  on  ground  which  we  will  call  estimation  profiles 
in  the  sequel.  Thus,  even  long  profiles  can  be  processed  very 
fast 

Fig.  4.3  shows  the  principle  of  the  moving  operator  R 
for  the  estimation  of  one  gravity  anomaly  in  each  step  using  two 

Flight  profilo  no  1 


0---0 — 


Flioht  DTofile  no  2 

C ~~  - -O^  - -C 3 


\ yy'v 
\ : 1 : ■ 

■ ' 1 1 / 

^ — / Interpolation  profile 


Fig.  4.3  Estimation  by  3 moving  operator 
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profiles  of  observations.  Obviously,  the  addition  of  nore  esti- 
mation or  more  observation  profiles  does  not  change  the  basic 
procedure.  We  can  estimate  point  as  v/ell  as  mean  gravity  anomalies 
by  this  method.  It  was  fc..nd,  however,  that  mean  anomalies  should 
not  be  determined  by  a mean  anomaly  covariance  function.  The 
smoothing  properties  of  this  function  will  cause  a considerable 
loss  of  information.  SLinkel  (1977)  has  derived  detailed  formulas 
to  estimate  the  loss  of  accuracy  for  a given  block  size  and  a 
given  covariance  function.  To  avoid  such  a loss  the  following 
method  has  been  adopted  to  determine  .mean  values.  Depending  on 
the  size  of  tne  block  3 to  5 estimation  profiles  are  used  at 
ground  level  and  gravity  anomalies  are  estimated  along  these  pro- 
files at  the  same  rate  as  taken  at  flight  level.  Thus,  in  our 
case  the  separation  of  the  point  gravity  anomalies  along  the 
estimation  profiles  is  2.6  km.  All  values  inside  a block  are 
averaged  to  obtain  the  mean  value.  In  this  way  the  accuracy  of  the 
point  estimation  is  maintained  and  we  obtain  an  averaging  pro- 
cedure which  is  linked  in  a simple  way  to  the  estimation  method. 

It  has  been  shown  in  Schwarz  (1976)  that  ag,  T , and 
T^  are  the  most  important  observations  when  estimating  gravity 
anomalies  from  east- west  profiles.  These  three  measurements  have 
therefore  been  simulated  in  each  data  point  at  flight  level.  To 
give  an  idea  of  the  accuracy  of  point  estimation,  fig.  J.A  shows 
about  3“^  of  an  estimation  profile  obtained  from  flignt  0''ofiles 
spaced  at  20'.  The  standard  error-  for-  the  yg-observations  is 
1 mgal,  that  of  the  gradiometer  observations  + 1 E . The 
normal  error  model  has  been  used.  The  heavy  line  shows  the  evact 
profile  of  the  simulated  field  at  ground  level  while  the  gravity 
anomalies  estimated  from  the  corrupted  data  are  represented  by 
dots.  The  agreement  is  very  good  with  a standard  error  of  about 
^ 3.3  mgal  for  the  point  estimation.  To  determine  how  much  of 
this  error  is  due  to  i nt erpo 1 a t i on  and  how  much  to  aownward 
continuation,  we  have  dete>^mined  an  interpolation  profil*^ 
the  same  data  but  this  time  located  directly  below  one  o'  ‘ne 
flight  profiles.  The  standard  error  reduces  to  jg  1.96  mga’.  "nu- 
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we  can  say  that  with  the  above  configuration  a measuring  error  of 
+ 1 mgal  at  flight  level  will  be  amplified  by  downward  continu- 
ation to  about  ^ 2 mgal  and  when  interpolated  between  two  oro- 
files  20'  apart  to  about  + 3 mgal  . A dense  profile  spacing 
will  therefore  improve  the  accuracy  considerably. 

Results  for  mean  anomalies  are  given  in  table  4.1,  The 
same  error  model  as  above  has  been  used.  The  number  of  mean  values 
estimated  in  each  case  is  given  and  although  the  samole  size  is 
small  in  one  case,  the  internal  consistency  of  the  different 
values  is  good.  For  comparison  results  obtained  in  the  corres- 
ponding accuracy  analysis  (Schwarz,  1975)  are  also  shown.  In 
general,  the  figures  agree  well.  But  the  accuracy  estimates  seem 
to  be  somewhat  smaller  for  block  sizes  below  3Q'x30'  and  larger 
for  block  sizes  above  30'x30'. 


Block 

size 

Prof i 1 e 
spacing 

Dumber 
of  values 

Simulation  study 
(mgal ) 

Accuracy  study 
(mgal) 

15'xl5‘ 

20' 

232 

+ 2.2 

+ 1.7 

30' x30' 

20' 

112 

+ 1.3 

+ 1.7 

l'"xl^ 

20' 

26 

+ 1.3 

+ 1.7 

l°xl° 

i 

1° 

78 

+ 4.6 

i 

+ 5,6 

1 I 

Table  4.1  Accuracy  for  different  block  sizes  as  obtained  f>'om 
simulation  and  from  accuracy  studies. 


The  influence  of  different  error  models  has  been  studied 
for  a number  of  cases.  Table  4.2  gives  results  for  mean  values  of 
15'xl5'  , a profile  spacing  of  20',  and  a standard  erroi'  of  the 

gradiometer  observations  of  m = + 1 E . 'lonmal  and  Markov 
models  have  been  comouted  for  different  variances  -.T  . The  '-'art-ov 
model  is  of  second-order  and  has  the  correlations  r,  = C.73 
and  r,  = - C.33  . 


Block 

Prof i 1 e 

Number 

Estimation  error  m_  in  mgal 

size 

spacing 

of  values 

(mgal  ) 

Norma  1 

Markov 

15 'xl5' 

20' 

232 

+ 1.0 

+ 2.2 

+ 2.3 

+ 2.2 

+ 2.8 

+ 3.7 

+ 5.0 

+ 4.5 

+ 10.7 

Table  4.2  Comparison  of  error  models 


The  results  are  remarkable  because  they  show  a signif- 
icant difference  in  the  general  behaviour  of  the  two  error  models 
once  the  variance  has  reached  a certain  size.  While  there  is  not 
much  difference  on  the  + 1-mgal  level  because  interpolation 
and  downward  continuation  contribute  the  largest  part  to  the  error 
budget,  the  difference  is  quite  obvious  on  the  + 5-mgal  level. 

For  the  normal  model  the  estimation  error  at  ground  is  smaller 
than  at  flight  level  while  for  the  Markov  model  it  is  more  than 
twice  the  size.  In  the  first  case  the  error  in  Ag  is  controlled 
by  the  accurate  gradiometer  data.  In  the  second  case  the  variance 
of  the  generating  process  is  enlarged  by  a factor  of  almost 
3 due  to  the  correlations  in  the  Markov  sequence,  see  formula 
(3.22).  This  raises  the  standard  error  of  = + 5 mgal  to  about 
= + 8.5  mgal  . The  actual  estimation  error  i 10.7  mgal 

is  even  larger  because  signal  and  noise  cannot  be  separated  as 
well  as  in  case  of  the  normal  model.  Thus,  results  will  deteri- 
orate considerably  if  the  errors  in  the  Ag-values  at  flight 
level  are  correlated  and  if  their  variance  is  not  extremely  small. 

The  situation  is  even  worse  if  the  data  are  corrupted 
by  systematic  errors.  Only  a few  examples  have  been  computed 
and  more  representative  studies  are  necessary.  But  a systematic 
influence  at  the  5-mgal  level  may  already  turn  the  estimation 
results  useless  for  geodetic  purposes. 

Finally,  some  investigations  have  been  made  on  the  in- 
fluence of  wrong  assumptions  in  the  covariance  function.  So  far, 
only  the  covar'iance  function  C(s)  with  the  essential  para- 
meters (4.<i)  has  been  used.  It  will  now  be  compared  to  the  co- 


A 


variance  function  C(s)  determined  from  the  simulated  data  with 
the  parameter  set  (4.5).  To  have  an  easy  distinction  the  latter 
will  be  labeled  'correct  covariance  function'.  Results  are  given 
in  table  4.3.  As  could  be  expected  the  correct  covariance  func- 


Block 

Profile 

Estimation  error  m^  in  mgal 

size 

spacing 

C(5) 

C(s) 

15'xl5' 

20  ' 

+ 2.2 

+ 2.1 

30'x30' 

+ 1.8 

+ 1.7 

•1  0 1 o 

1 X 1 

+ 1.3 

+ 1.0 

rx  1° 

+ 4.6 

. 

+ 6.2 

Table  4,3  Comparison  of  different  covariance 
f uncti ons  . 

tion  ^(s)  gives  slightly  better  results  for  all  block  sizes  when 
a profile  spacing  of  20'  is  used.  The  differences  are  not  large, 
however,  and  we  can  conclude  that  the  choice  of  the  covariance 
function  does  not  affect  the  estimation  very  much  as  long  as  the 
essential  parameters  are  reasonably  close  to  the  empirical  values. 
This  confirms  results  published  in  Moritz  (1976).  It  should  be 
kept  in  mind,  however,  that  these  results  are  valid  for  isotropic 
covariance  functions  only  and  may  not  hold  in  mor'e  difficult 
cases  . 

The  last  line  of  table  4.3  which  shows  the  mean 

*or  a pr'ofile  spacing  of  apparently  contradicts  the  above 
conclusion.  In  this  case  the  C(s)-function  gives  much  better  '-e- 
sults.  The  most  probable  explanation  is  that  the  data  are  not 
dense  enough  to  get  a useful  estimate  from  the  ir(s)-function.  The 
correlation  length  T is  smaller  than  orie  half  of  the  profile 
spacing  and  it  has  been  found  that  results  may  become  very  poo*' 
in  such  a case.  Thus,  estimation  with  a non  o o t i m a 1 c o v a r i a n c e 
■unction  may  in  some  cases  give  better  results. 


5.  Brief  Description  of  Programs. 

The  programs  and  a sample  computation  have  been  added  to 
the  report.  Since  detailed  comments  are  given  in  the  subroutines 
we  will  present  only  a brief  survey  of  the  programs  and  point 
out  a few  restrictions  in  the  present  version. 

Basically,  the  program  computes  point  or  mean  anomalies 
at  ground  level  from  measurements  at  flight  level.  Terrestrial 
and  satellite  observations  may  also  be  added.  It  has  been  assumed 
that  some  preprocessing  of  the  airborne  data  has  taken  place,  and 
that  besides  5 second-order  gradients  the  value  of  Ag  is  avail- 
able in  each  observation  point.  In  the  present  form  a total  of 
7 observations  per  point  can  be  read  but  fewer  are  possible.  All 
data  should  be  available  on  a direct  access  file.  The  subroutine 
ZINF  selects  the  observation  profiles  needed  for  the  estimation 
of  a specific  point.  The  subroutine  DATS  reads  these  profiles  and 
deletes  those  measurements  which  are  not  needed.  In  this  way  the 
core  storage  requirements  can  be  kept  relatively  small.  Further- 
more, the  simulated  data  from  the  file  are  corrupted  by  one  of 
the  error  models  in  DATS.  This  part  of  the  program  must  be  deleted 
when  handling  real  data. 

The  number  of  estimation  profiles  for  a specific  mean 
value  is  fixed  in  VMEAN.  We  have  used  5 estimation  profiles  for 
blocks  larger  than  40'x40',  3 profiles  for  blocks  larger  than 
5‘x5'  and  smaller  than  40'x40',  and  1 profile  for  blocks  of 
5'x5'  and  point  values.  If  more  than  5 estimation  profiles  are 
used  the  correspond!  ng  D IME'IS  I OW- s ta  temen  t s have  to  be  changed. 

The  covariance  matrices  are  set  up  in  the  subroutine  PLAY  which 
uses  the  subroutine  COVAX  to  determine  the  individual  covariances. 
The  estimation  according  to  formula  (4.1)  is  done  by  the  subrou- 
tines FIX  and  SCAN.  The  subroutine  COMPA  compares  these  estimated 
values  to  the  true  values  of  the  simulation  which  are  again 
stored  on  a direct  access  file.  The  error  computations  are  also 
performed  by  this  subroutine.  Again,  this  part  of  the  program 
must  be  deleted  when  handling  real  data. 

Three  of  the  subroutines  used  are  not  listed.  The  sub- 
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routine  COYAX  has  been  published  in  (Tscherning,  1976).  Three  cor- 
rections which  have  been  communicated  by  the  author  are  listed. 

The  subroutines  DSIflV  and  GAUSS  are  available  in  the  Fortran 
Scientific  Subroutine  Package  which  can  be  found  in  IB’-l  publi- 
cations. The  first  one  inverts  a double  precision  matrix  stored 
in  a one-dimensional  array.  The  second  one  generates  normal  de- 
vi ates  . 

In  its  present  form  the  program  accepts  flight  profiles 
in  east- west  direction  with  4200  observations  each.  This  corres- 
ponds to  a profile  length  of  lo”  if  the  flight  speed  is  500  knots 
and  6 measurements  are  taken  at  each  point.  From  a total  of  11 
profiles  a maximum  of  5 profiles  is  selected  for  each  estimation. 
Since  the  storage  requirements  are  strongly  dependent  on  the  length 
of  the  profiles,  a subdivision  of  longer  profiles  might  be  con- 
sidered. Because  of  the  use  of  the  direct  access  file  an  increase 
in  the  number  of  profiles  is  not  critical  as  long  as  the  maximum 

number  of  profiles  in  the  core  storage  is  kept  to  5. 

The  speed  of  the  data  processing  depends  mainly  on  the 

number  of  flight  profiles  used  in  each  estimation  and  only  to  a 

smaller  extent  on  the  block  size  of  the  mean  value.  As  an  examole 
mean  values  covering  an  area  of  5'^xl3'^  have  been  comouted  using 
the  data  configuration  and  the  covariance  function  given  in 
section  4.  With  an  average  of  3 flight  profiles  for  each  estimation 
about  30  seconds  of  CPU  time  are  needed  to  determine  75  blocks 
of  I'^xl'^  ; for  5 profiles  the  time  required  increases  to  about 
50  seconds.  If  we  estimate  300  blocks  of  15'xl5'  instead,  another 
10  T must  be  added  to  the  time  estimates.  Thus,  even  large  amounts 
of  data  can  be  processed  in  a relatively  short  time. 
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6.  Conclusions 

The  simulation  studies  presented  in  this  report  show  that 
1 ea s t- squ a res  collocation  offers  an  adequate  model  to  estimate 
gravity  anomalies  from  airborne  gradiometer  measurements.  The  pro- 
cedure is  simple  numerically  and  allows  to  handle  large  amounts 
of  data  with  regular  requirements  on  core  storage  and  small  demands 
on  computer  time. 

The  deviations  of  the  estimated  gravity  anomalies  from 
their  true  values  agree  well  with  the  estimates  obtained  from 
correspond i ng  accuracy  studies  (Schwarz,  1976).  It  should  be  noted, 
however,  that  correlated  errors  in  the  measurements  will  strongly 
influence  the  accuracy  of  the  results.  A second-order  Markov  sequence 
has  been  used  to  model  the  error  process  along  the  profiles.  De- 
pending on  the  size  of  the  correlations  and  the  variance  of  the 
process,  the  mean-square  errors  will  more  than  double  as  compared 
to  the  uncorrelated  case.  Similarly,  a bias  in  the  data  will  im- 
pair the  accuracy  of  the  results  considerably. 

The  simulation  of  gravimetric  quantities  from  a point  mass 
model  is  considered  in  the  spectral  domain  and  conclusions  are 
drawn  with  respect  to  the  resulting  fields.  In  order  to  represent 
adequately  regional  variations  of  the  gravity  field  as  well  as  the 
local  behaviour  of  the  second-order  gradients  the  medium  and  the 
high  frequency  part  of  the  spectrum  must  be  modeled  equally  well. 

This  can  only  be  achieved  by  using  several  planes  of  generating 
masses  at  different  depths.  In  many  cases  a model  with  two  planes 
may  already  be  sufficient. 
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Appendix  A:  Computer  Programs 


SIMULATIO.N  SruOY  AIR9GRNE  GRAOIO''ETRY 

The  PROGRA''  computes  point  and  mean  gravity  Ar.nf-iLJcS  ON  GROJNO  FRr)'< 
aerial  gravity  a, no  O-'ACIOsETRY  •‘EaSG''.l“EMS  ago  FROR  GEOIG  iNECR’-'AT  ion 
PROviDEu  oY  Satellite  Ai_ri«cTRY.  t-e  conputeo  valjES  are  compareo  tq 
•TRgE'  VAl'JES  GENcRaT.'O  “.V  A MASS  anomaly  "JOEL.  DEVIATIONS  ARE  GIVEN 
FOR  ALL  CGMPJIED  VAlOLG  JND  TmE  EFFECT  JF  DATA  FE.RT  liR  b\ T I JMS  VARIOUS 
xiNOS  Can  3E  uEie-^mineo.  at  The  present  stage  the  Program  is  rfst^icteo 
TO  ErST-RcST  PruFILES 

implicit  real  -'3 ( A-M . ^-Y  ( , LCGIGAHL ) 

COMMON  /TICn/  3P(12)..:?llE)»0P(l?i.a(lb),Cll6).ullS).FR{l!j»,FCIlH) 
UTRaNS  ( 12)  .SR  I S ) ,EO,;CO.  V’/M.SMI  .Ran  ,z  ( 5,  JEOO  I . YMR  I 12)  t<v<p  I 1C  ) . lOP 
2I13.r),R<(lo),RR.Mil6).[D(16.7)  ,NN,MM.N“.N.NN,  IM  I 1 E 1 . NNV  , NV  . M V , I C<  . 
3rjj,.Mvv.M''*v  (IOI.itz.ifi 

COMMON  /CMC  JV/CI  (l?),v.Rl5U.SIG«AO(  300  I . SIGMAI  300)  . )^I  ( 25  ) ,N1  .LOCAL 
common  /ShT/  C2( 3, S) 

0IM£))SI0N  AA  I •.o,  Ao  ) , 93  I 23.  "6  ) , P.B2  123.23)  .G0(  10)  .GCO  i 10)  . GN(  ID  ) 

DI"£NSIJN  AAA(  UOO)  ,RR(e,12.6) 

FOMl valence  ) AA(  1, li  ,A4a(  I ) ) 

C<-  Ta  &M,  RE/3.  G 901  m.6  371.003/ 

*.DO.Ol.D2/'O.OuO.  1.003,  J.OGO/.F  1/2  .lAl  5926  53500/ 

•RE  MEAN  RAO  US  OF  THE  EARTH 

GM  RR30UCT  Or  THE  GRAVITATIONAL  CONSTANT  AND  TnE  MASS  OF  THE  EA’TH 
the  common  area  /tier/  is  Used  to  TRANSFER  DATA  FROM  THE  -AIN 

PkOGPAm  to  all  other  SUjP.C'JTINES.  The  common  area  /C'COV/  is  used  'j 
transfer  GATA  FkOM  T-E  MAiN  PROGRAM  VIA  THE  SUBROUTINE  PLAY  TO  THE 
SUHRCuTINE  C'VAX. 

CO  = 1. 030/57. 2V5779E1D0 
COO  = CO/60. 

MA  = 66 
MB  = MA/2 
00  110  1=1. MA 
00  110  J=1.MA 
00  110  K = 1.,M3 
AA(I . J ) = 0. 

110  Eu(n.I)  = 0. 

DO  111  I =1.0 12 
RRII.3)  = IOC. 

00  111  0=1.2 
Rkl  I . J ) = 0. 

DO  111  R=H.o 

111  RR( I .<  ) = 0. 


INPUT  COVARIA.Cb  F'JnCTIUN 


100  READ  (5. A)  S.a.yK 5) .K’.<3,N.L0CAl 
9 FORMAT  ( 2U1h.  7 . Al  E.L.;  ; 

IF  (N.NE.O)  GU  TO  101 
LOCAL  = .TRUE. 

N = 2 

101  N1  = '..I 

iFI.NoT.LuCALlREADIN.lOjlSIGMAO:!).  I 


1.  N1  ) 


I 3 FQRyAT  I i3R‘>,2l 

IF  1 .NOT  ,l2C  il.  ) ViKl  Te  ( if  7 ! I t IG-AC  t : ) , ! = If  Sit 
7 Fo^.SAll*  f’r'lK.UAL  FNC-aLY  Dc'irl  £ E -V  A R I A n C c S IS  oSITS  OF  -OAlASS:*! 
Vf25(12Fi.;/)  I 

c 

c s 

c 

C A 

C Klt5) 

c 
c 
c 
c 

C K2,<3 

C N 

C LGOAl 

c 
c 
c 

C INPUT  FLIGHT  PROFILES  Ar.O  SEASURi NT  S 

C 

READ  ISfll  N-lf  VHV.fVI  fA:<L 

1 format  1 I5,20l2.2f 15 ) 

00  2 I=I.SN 

READ  (5f?)  o(  1 1 fC  I : ; f 01 1 1 f GNt  1 1 f sot  I ) f Gool : ) ,k:<  1 1 ) , rrk.!  n , 
ii  I0( If ji , j = i t n 

3 FOS-AT  (607.D,X5,JI3) 

2 Oil)  = DIDFIOOO. 

C 

C NS 

C RK 

C KKK 

C AKL 

C V M,1 

C V I 

C b 

C C 

C 0 

C 10 

c 
c 
c 

C 2 

C 

c 
c 
c 

C ON 

C GO 

C GOO 

C 

C INPUT  iNTcPPOLAT lOS  PkOFIlES  ASO  MEaS  ANOMALY  9L0l<S 

r 

READ  i5.5(  loi 
5 format  (a;5.2P!0--> 

SNS  = N 4*  I 
SM  = 

K£  = 0 

00  2 7 : =NNSf 'JX 

READ  'Stilt  2(  I t .:t  1 t f'M  It  fASt  n 

II  Fl,R-aT  (307.2fI3l 


^ 


NUMOER  of  FLIGHT  PROFILES  CftSERVEO  | 

NU'AOER  OF  POINTS  =>ER  PROFILE 
NUM3ER  or  ■'FaSL'REHc'.TS  per  PjIST 

NUMOlR  of  REOOkOS  to  5i  READ  FaO“  EACH  PLIGHT  FILE 

MEAS  VELOCITY  OF  AIROP.AFT  IS  RNOTS  PER  hJ.jr 

MEAN  interval  oET.ffc-=.  05  SS  R V A T 1 OftS  IN  SECOSOS  0=  Tr'j 

LATITUDE  OF  ISITIAl  JOIST 

LOSGITuOE  JF  I.-flTIAL  °JINT 

ALTITUDE  OF  INITIAL  P I S T 

SUMRtPS  specifying  '•’EASJRFhCnTS  AS  Is  TSC-'£RnING  (17  7 5) 
possible  NCHoERS:  ItjfStiO  ’'j  Ia 

SEOJENCE  OF  ,U-'3cRS  A S C ;■  S')  I NO  . I = LESS  Than  y - E a sjr  £ vE  N T S t 
THE  RFmaIMNO  F'tSITMNS  -JST  9k  FILLED  5Y  O7ROS. 

MATRIX  storing  T-iE  0.j  5E  R VA  T I OSS  . c JL  h SOSTaINS  jSc 

PROP  I LE.hEaSURE"ENT5  are  GT'JREj  fl'iST  s'  pjInT  Ii  'r-s 
SEUUENLE  OESCRIOEO  ?Y  IO.p-OInTS  aRc  s.’uEfEO  aS.'RGISS 

TO  asccNoing  lung; tude, profiles  accoroino  to  7es::.^i.o 
latitude. 

Variance  of  gsuijal  uscuLirijss  is  heter  '■ouareo 
VARIaSlE  Or  DElTa  G IN  'GAL  SGgaRFO 

VARIa'CE  Or  SECO.no  OE  .1  IV  AT  I v E S IN  EU  SJUAPEO  | 


square  of  T-HE  ratio  3ETGEP,N  radius  of  The  aOER-iAHfcP  sphere  aS^ 
mean  RADIUS  0=  The  Earth 

COfSTA.-.T  paCTCk  of  GRAVIT’  A.NOMAlY  DEGREE- v a k I a r;C  E ••■□Ccl 
SUY3ER  OF  C£..RPE-/ARi  ANOE  MODEL  AS  SP'CIFlbO  IN: 

TSCHER.NlsG)  IvTil  :C^'vaR1a.iCc  EXPRESSIONS  FOR  SrCGtO  AND 
LOi.ER  ORjER  o=^:.Arrv£S  J=  The  ano.malous  pctcniia. 

OSJ  .REpnp.T  NO  225 
P0SS:3L=  S.'.Er.S:  1,2.3 

I.NTEGErl  n2  .-.no  Sj  IS  '-GRPOlA  (17)  IN  T SC  hERN  ING  ( 1 R 7 5 ) . ; 

DEGREE  VARIANCE!  J=  TO  ASO  InCLUCISG  DEGREE  N ARE  -ITh£R  SET  -■  'j 

ECJAL  to  OERO  or  JsPLACEO  s'  PHPIP.ICAL  OtO'EE  V ah  ; .•,-.CP  G . IS  T '£  A ' 

FIRST  Case  T,';P  LO.SICal  va’IaELs  lOCAl  MUST  ^E  .'aUc.  In  IHc  ; j 

seco.so  CASE  .paLSE.  THE  l'PI.RIOAL  SRAVITY  ASL-‘ALY  JEGREE  i 

-VARIANCES  “UST  EC  R£AU  IN  J.ITS  OF  MGAL  S.jAREO.  j 


r'  n o oooononononon 
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27  <e  = KE*<K(I) 

IT2  = ITZ-'f 

READ  CIFI,32)  ( Z (1 , J I t J = 1 1 1 TZ I 

32  format  lAAAl  --i 

X = 0 t 

QO  31  1 = 1,  ITZ.A 
K = K*  1 

RR(K,A I = Z( 1 . I ) 

RR(X,5I  = Z11,I*1I 
31  RklX.ol  = Z1 1,1*2) 

IFI  = ITZ/‘* 

lOl  -1  POINT  VALUE 
1 MEAN  values 

MM  number  of  interpolation  profiles  - POINT  OR  rt'^AN  VALUES 

kEOUIREO  NU'*oER  of  INTERPOLATION  POINTS  PER  PROFILE  - POINT  >i 

OR  mean  values  • 

IFI  file  NUMjER  where  “EAN  values  of  TRUF  ''GDEL  are  STORED 

ITZ  number  Or  MEAN  VALUES  STORED  IN  FILE  MFl* 

VM  SIZE  OF  MEAN  ANOMALY  SLOCX  IN  MINUTES  OF  ARC.  MUST  BE  ZERO 
FOR  P0I‘!T  values. 

SMI  MEAN  profile  SPACING  IN  MINUTE.V  OF  ARC 

ALL  OTrtFR  gUANTITIES  AkE  AS  OEFINEO  FOR  THE  FLIGHT  PROFILES 

IF  ( lOI I 2A,22,23 

22  WRITE  (6,25)  101  « 

25  format  (///2x, 'Parameter  ioi  cannot  be  zero. execution  interrupted' 

1,/// ) 

CO  TO  301 
2N  VM  = 0. 

DATA  CONTROL 

23  WRITE  (6.1A) 

lA  FORMAT  (/SOX.'O  ATA  CDNTRQ  L',//50X) 

WRITE  (0.12)  S.A.M ( 5 I ,X’,<j,N,lOCAL 

12  FORMAT  ('  parameters  SPECIFYING  ThE  OEGR FC- V AR I aNCE  “uOE L ' . / / , I OX . 

I'S  =',F1A.6,/10X,' A = • .F13.  A,/ ICX  ,•^^  ( 5)  = ' , I / I OX  , ' K2  =',I7  ,/10 
2X,'X3  ='.17  ,/lQX,'N  =',I3  ,/ 1 OX , 'LOCAL  =',lA/| 

WRITE  Io.2C)  VMM.V I , VM. VM.SMI 

20  FORMAT  (//ZX.'mEAV  VELOCITY  OF  AIRCRAFT  IN  KNOTS  PFR  H 0 j R : ' , 9 . 2 . / 

12X,'“EAN  INTERVAL  BETWEEN  MEASUaING  POINTS  IN  SECONDS  OF  TlMEt'.F? 

2.2,/2a,  ' SI  ZE  OF  «EAN  ANu“AL  T bLuC  K:  • . F 5 . 0 . ' * • , 5 .0  , • “ I UT  S ' , / 2X  , 

3'MEAN  PwTFILE  SPaC ING: ' , F7.2, • MINUTES'.///) 

DO  I?  r=l,NN 
WRITE  (6,18)  I 

13  FORMAT  (2X, 'FLIGHT  P»0=tLE  NO ' , I 3 ) 

WRITE  (6.15)  b( I ) ,C( I ) ,0( I ) ,KK( I) 

15  format  (2x, 'Initial  point  of  profile  ^ atituof: • . = 7.2. ' lCngitu 

IDE:',FY.2,'  altitude:' ,F7.0,/2X,  'number  oF  POINTS  IN  PROFILE:',!- 
2) 

19  write  (6,16)  GN(  I ) ,C0(  I ) ,GDD(  I ) , ( n ( I , J i ,0  = 1 ,7 ) 

16  FOPMAT  ( 2X, 'VARI ANCE  CF  GEOIO  UNC JL A T I ON S : ' , F 7 . ? , ' '=  DELTA  G:' 

I, FT. 2,'  OF  SELCNO  ORDER  DcR I V A T I V £ S : ' , F 7 . 2 , / 2 X , ’ S PF C I E I - A T I ON  No 
2M0ERS: ',713/) 

WRITE  (6,951) 

WRITE  (6.951) 

951  FORMAT  (2A/) 

00  A 1=NNN,N" 

Ik  = I-NNN*l 
WRITE  (6,17)  IR 

17  FORMAT  ( 2;;  , ' INTERFSlAT  ICN  PROFILE  NO', 13) 


0(11  = D( I )-100Q. 

<►  WRITE  (6tlSI  S ( I 1 .C(  t I .C(  n tKK(  I ) 

computations 

WRITE  16, 9 ‘5  1) 

WRITE  (0,9511 
WA ITE  {o,',5£) 

952  format  (/SOx.'C  OmouTATION  S*,//) 

covariance  model 

RB2  = R£=RE^S 

KI(j)  = <2 

M ( 5 ) = <3 

LMOOEL  = N.cI.O 

CK8)  = A9RS2*  I.GD-IO 

Cl ( 10)  = S 

CALL  COVAX 

DcTtR.MINAT  ION  OF  FLIGHT  “ROFILES  NECESSARY  FOR  INTERPOL  A TI  0 , 
NUMbEKEO  ICnE 

CO  30  1=1, NM 
3(  . ) = 3(1  )*C0 
30  CO  = C(  I )RCJ 

WM  = V)'.M5VI/360Q. 

CALL  VMEAn  (VM) 

IF  IV.M)  3J1,2S,23 
23  ICXE  - NNN 
VM  = V^RCOO 
VVM  = VyMRCjO 
IMM  = 1 

300  CALL  2INR  (IC.nE.VM! 

IF  ( IT2)  302,29,2, 

29  IC  = ICXE-NN 

WRITE  (6,9531  IC  , ( TR ANS ( J) , 3= 1 ,NJ I 

953  format  ( //2X. • interpolation  PROFILE  NO*, 13,*  OSES  FLIGHT  PROFILES 
INOS* ,5F5.0,// I 

COMPUTATION  Or  aUT OC 0 V AR I ANC E “ATRIX  OBSERVATIONS  CXX 


303  CALL  PLAY  (AA,3B,-I,MA,M3) 

“U  = 0 

00  23o  1=1, NO 
230  Mu  = MJ»<.<P(  I )*<.X-'.?l  I 1 
JLL  = NJ»MV 
JJ  = 0 

OO  232  I = 1 ,NJ 
XRR.s  = XF?i:) 

00  232  ir=l,<vYK 
00  232  „ = 1 , 7 

IF  ( ILP(  I, J! -1  I 23  3,23a,235 
233  CONTINUE 
GO  70  232 
23<,  JJ  = JJ*1 

IN  = TRA'IS(  I ) 

4A(JJ,JJ)  = A A ( JJ, Jj ) »ONI I N) 
GO  TO  232 

235  IF  (109, [,21-31  233,323, 32L 
323  JJ  = ^J«i 

IN  = T'ANS(I) 
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AA(JJiJJ)  = AAIJJ. JJ ) ♦GO! INI 
GO  TO  232 
32A  JJ  = JJ*i 

IN  = T«ANS( I ) 

AA(JJtJJ)  = Aa( J J, JJ ) *COO( INI 
232  CONTINJt 

INVERSION  Of  CXX 

CALL  ORDER  1 A A , MA , NO , - 1 1 
SOF  = I.E-50 

CALL  OSINV  ( AAAf  >'U.iOF»  lERI 
IF  (lERI  95StN5e.,'yf5 
956  WRITS  (6,9731 

973  format  ( //2X, • INVCRS ion  OF  COVARIANCE  MATRIX  OKAY", //I 
GO  TO  2A3 
955  WRITS  (5,9791 

979  FORMAT  ( // 2X , • I NST A3  I L I T I E S OCURREO  DURING  THE  EXECUTION  OF  TmS  IN 
IVERSIDN  RROGRaM.  EXECUTION  I NTEKR U’ TE J . • , // I 
GO  TO  302 

293  Call  order  (aa,ma,mu,ii 
IF  (JJ-MJI  230,237.236 

236  WRITE  (6,299)  JJ.M’J 

299  FORMAT  (/20X, 'NUMSER  0=  ROWS  IN  CXX-MATkIX  ERRONEOUS  JJ = • , I 3 , • “U= • 
1,13/1 

237  00  231  I=l»MJ 
DO  231  J=1,M0 

231  AA( J,I I = AA( I,JI 
MUM  = MV 

COMPUTATION  OF  PREDICTION  MATRIX 

CALL  play  (AA,3e,0,MA,MLl 
DO  2 39  J = l,M'J 
DO  290  1=1, HUM 
SUM  = 0. 

DO  290  JJ  = I, MU 

SUM  = SUM» 33( I , JJS »AAl JJ , J 1 

290  BS3( 1,1)  = SUM 
DO  291  11=1, MUM 

291  AA( I I , J 1 = E33(  11,11 
239  CONTINUE 

COMPUTATION  OF  SIGNAL  (POINT  OR  mean  VALUES) 

I2I  = -i 

Call  DATS  (KXL, S.N, GO, ijO 0,121) 

CALL  SCAN  (AA,Ma,;“M,VM,RR. ICAE) 

DO  292  I=I,MJM 
P'O  2 92  J = l,MJM 
SUM  = 0. 

CO  250  JJ=1,MJ 

250  SUM  = SUM, aA(  I . J J I P3J( J, JJ ) 

29  2 BQ3(I,J)  = ^U.M 

COMPUTATION  OF  COVARIANlE  MAIRIx  UF  SIGNAL  CSS 


CALL  PLAY  (AA,33,l,MA,Mj| 
DO  307  I = I..'U.-I 
DO  307  J=l,MUM 
307  AAl J,1 I = AAI I,JI 
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C COMPUTATION  Of  ERROR  COVARIANCE  «ATRIx  jF  5I3.NAL  ESS 

C 

00  2R3  1=1. MUM 
DO  2 A3  J = 1.NUM 

243  333(1,0)  = AA( I,J)-S33(I.JI 
ICiAE  = ICRE*l 
IF  (ICKE-.VM)  300,300,302 
302  ICRE  = NNN 
C 

c print  results 
c 

00  306  1 = 1, <E 
306  RR(  I , I ) = RR ) I , 1 )/CJ 

00  303  I=‘.NN,NM 

308  3(11  = 3(  n/CU 
SUN  = 0. 

DC  310  1=1, NUN 
313  Sun  = SuN.BuO ( I,  I I 

IF  (NjN-U  311,311,312 

311  Su-S  = 1 
CU  TO  313 

312  SUNS  = njm=( 

313  SNI  = DSl.lT  (SUN/SUNS) 

CALL  C3.N=>A  ( Rr.rE,  Sj“',  SUNS  : 

V.N  = VN/CCO 

IP  ( IJI  ) =56,  =69, OS'? 

969  WRITE  (6,d)  V.N,V« 

6 FQRnAT  (//4^X, 'MEAN  values  FOR  dLQCKS  OF',=5.0,'  »',=5.0,'  NISUTES 

I’,//) 

(SO  TO  970 
968  WRITE  (6,972) 

972  FORMAT  ( // SO  X , • PC  I NT  VALUcS’,//) 

970  WRITE  (6,?71) 

971  FORMAT  (luX,*  NO  • , S X , • L AT  I TUOE  • , 5X  , • L CNO  I ’’UOE  ' , S X , • C . w PUT  E 3 • , ? X . • 
1TKu:-*,7X,'0IFFE.'wE',C£S'.  sX.'pERCENTaSc'./mTX.’VALjES'.Tx.'VALjES', 
210X,  *>,-1  • , 7x, 'OF  necessary  NEAS.*,/) 

1 = 1 
JL  = 2 

jLL  = <r.  ( ISXE  ) 

304  WRIIE  (6,6-:)  I , 3 ( I C<E  ) ,RR(  I , 1 ) ,RR  ( I , 2)  ,RR  ( : , S ) ,RR  ( I ,s)  , RRi  : , 3 ) 

68  format  I r 1 3, if  1 a.2 ,oX ,F7.2.9X , F7. 2, 7X, F7 ,2,r Is.O ) 

30  66  1=Jl,^L,. 

66  WRITE  (£.,67)  I , ( I , 1 ) ,RR  ( : , 2 ) ,RR(  I , 6 ) ,RR(  I .6  ) ,.RR  ( I . :•  ) 

67  d^6;,AT  (IL3,  F 26.2 ,6X  ,'-'7.2 ,9.<  , F7. 2,  7X,  FT.2.F  16.2  ) 

ICXE  = ICnE.I 

JL  = JLL*2 
I = I.l 

JLL  = JLL.<\(  IC.vE) 

IF  (12<E-NW)  3Cm,33u,3'3E 

305  WRITE  (6.339  ) ' u “* , < E , SUM  S , SN  I 

309  =ORnAT  (///:.  y ,•  SU.N  SF  j are  S ; • ,F  U . 2 , /2X  , • S TanCARO  DEVIATION  =kVa 

l',I5.'  nEaS^RcMENTS:  ' , = 1 3.  3 , • ‘'o  A L ' . / 2 x , ‘ F ')■,  m a i.  sta.Sarj  O'^vIa’ISn 

2 FROM  COLL  3S  A T I CN  : ' , F / , 2 , ' “SAL  ' , //' ) 

3CI  STOP 
END 


SUMRCUYIN;  vNcAN  ( V“ ) 


C 

C 


DETERN 


UANTIlriS  FJR  Trie  “EA.'w  VALUE  C.'NPuTaTISN 


r»oo  nooonr»  n nnnononnoo 
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VM  MEAN  SIDE  length  TF  DLOCK  IN  MINUTES  OF  A^C 

RAN  ZONE  Of  INFLUENCE  ON  EaCFi  SIGE  Of  I NTcRPOL AT  I ON  PROFILE. 

N.N7  NUM3ER  CF  POINTS  rCR  ONE  MtAN  VALUE  AND  UNc  FLIGHT  PROFILE 

MMV  numoer  of  points  On  each  flight  profile  used  to  SFT  up  The 

ESTIMAIIJN  operator.  mhV  must  EE  ODJ  ANO  GPEaTER  ^han  I. 

MW  separation  of  the  aoqve  points  in  units  GF  VV“ 

MV  NUMBER  Or  IMEkPOLATIQN  PROFILES  Fu»  “EAN  VALUE  OF  SIZE  VI 

NV  NUMBER  Or  POINTS  O.N  EACH  I NTE  P POL  AT  I ON  PROFILE  USEC  FOR  THE 

ESTIMATION 

IMPLICIT  REAL  Po ( A-H , 0-Y I . LCG I C AL ( L I 

CQM.NOn  /UC.</  BPI  12  I .C®!  UI  . op  1 12  .51  IS)  tCl  IS)  tO(  if  t.FE(  Is  I ,FC  ( IS) 
lfTRANS(12I.SR(5),CO,CCC,VVM,SMI,RAN,Z(S,A200'.RRP(l2I.Ki.K;P(l2l.IOP 
2(l2r7),^Y(la).<^‘'.(l^),I01lS,7)  .NN  , M.M  . N.M,  NNN.  I N I ( I 2 I . N.NV  . o . MV  . I CM  . 

3N J fMVVfM^VtlOIflTZtlF) 

MMV  = 3 
MW  = 5 
NV  p 1 

RAN  = I.59VMPCOO 

SMI  = l.5*SMI*COQ 

IF  IRAN-SMIi  901. VO.;, 902 

901  RAN  = SMI 

902  SMI  = 1. 

IF  ‘VM-AO.)  910,910,311 

911  MV  p S 
ICR  p -I 
GO  TO  912 

910  IF  lVM-5.)  919,913,913 
913  MV  p 3 
ICR  = 0 
GO  TO  312 
91A  MV  = I 
ICR  = I 

IF  (VM)  912,915,912 
915  V.M  p VVM 
NNV  p 1 

912  RETURN 
ENO 


SUBROUTINE  ZINF  (ICkE.VMI 

determines  jHp  profiles  to  oF  used  FOk  IHE  I NTFRPOL AT  I ON  OF  A SPECIFIC 
POINT  usi.s  The  latitude  ,-,lat  of  This  point  a iO  Th*  Oua'.titv  ran  as  f;teo 
BY  the  SJuRwUTINE  VMEAN.'hE  PROFILE  NuMoESS  ARE  STOkSU  IN  ThF  ARRaY  TRaNS. 
ICRE  number  of  interpolation  PROFILt 

implicit  real  *8 1 A-H,o-r ) .logical (L  I 

COMMON  /UZ</  3r(i2l,oP(UI,0P(12),c(lS».C(ISI,j(lS).F3(!sl.F:(Iol 
l,TRANS(l2I.SR^5|,C3.CCO,W“,SMI,RAN,Zl5.•.^0J|,K^^(l.;),5^kP(i^|,iJP 
2(12, 71, R<lls),Rr.R(l6l, lulls, 71  ,NN  ,'*M,  NM,  N)iN,  I Ni  ( I 2 I ,NNV  , Nv  , “V  , ICR  , 

3NU,MW  ,MMV,  lUI  , I T2  , IF  ! 


SELECTION  OF  FLIGhT  PROFIlES  FOR  IN  TF  R POL  A T I CN  PRO=IL£  N'JMaERtO  IC«t 
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NNV 

= VMP3C 

0S( 3LAT )/VVM 

IF 

1 HVV-5 1 

419,419,420 

419 

MW 

= 5 

420 

IF 

1 NNV-1  ) 

421.421.422 

421 

NNV 

* K.S(  IC.<E) 

422 

ITZ 

= 1 

STI  = 1VV 
STT  - STI»VV'* 

STI  = 

00  ‘»10  I-l  .10 

410  TRA.-iS(I)  = 0. 

BLATI  = 3LAT*sAN 
SLAT2  = 3LAT-RAN 

1 = 1 
J = 1 

401  3L=  BID 

IF  (3L. or. BLAU)  00  Tj  00 
IF  t3U.LT.cl.AT2J  OC  ^ L,  400 
TFAr.StJ)  = I 
BPU  ) = Bin 

CPtJI  = 3l0-STIpSTT/0C0SI3( 1 n 
OPtj)  = Oil) 

MUPIJ)  = 

lUi^AP  I J ) - <K.<  1 I I 

FStJ)  = 0. 

FCl J)  = $TT/OOOSICt I I ) 

DO  409  A=l.7 
409  lOPtJ.K)  = lOlI.A) 

NJ  = J 
J = J»\ 

400  I = 1*1 

IF  II-NN)  401,401.402 

402  IF  INJI  413, 4,3.412 

413  WPITF  la.4141  BOAT 

414  FO=“AT  (///:.X,’NQ  Flight  P9C=IL£S  in  Tnt  OBFINBO  PE  G I ON.  1 4 t P ? ?L  A 
IIGN  PPOFIuf  «ITrl  Tit  LATITCuc  3l  A T - • , r 7 . 2 . ' LA.nNQT  EE  OtTEl-INcO 
2./// ) 

IHH  = I«H»<«.P  I ICKc  1 
IC4E  = IC«,E»l 
IF  lU4c'‘(*M  416. 415,417 
417  IT2  = -1 
GO  TO  413 
C 

C PATTEP't  OF  PUePPOLATT  ON  POINTS  =09  SPECIFIC  4EaN  VALlE 

C values  used  fop  I-.TcPPOlAT  ion  of  one  PPjFILE  ape  ST_TcD  In  /BSO/ 

c 

412  STI  = HV 

ST  = V-/STI 
J = SJ.l 

IF  lies)  403*404, >.05 

403  3?(J)  = 3LAT.2.SST 
ePlJ*U  = lu.  AfST 
?P(J*2)  = '.'wAT 
BP(j.3l  = 3lAT-ST 
SP(J*4I  = iLAT-2.;ST 
GO  TO  406 

4C«  ?P(J)  = ^.AT.ST 
?.P'J*l)  = 3wlT 
eP(J»2l  = 3LAT-ST 
GO  TO  436 
4J5  2P(jl  5 BLAT 
40*  00  407 


I 


Ik 
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J = NJ«K 
CP(JI  = ILO 
OP(Jl  = Hhl 
KKPlJl  = 

PBUI  i 0. 

FC(J)  = STT/3COS(bLAT) 
KKKP(J)  = I 
lOP(Jtl)  = 3 
00  407 

407  lOPlJ.UI  = J 
418  RETJR^4 
END 


I 

I 

I 


I 

1 


C 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


SUBPQjirJE  PLAY  (AA,BB.MQL0tYAfM9l 
COMPUTES  covariance  MATRICES 

MOLD  -1  AUTOCOvARIANCE  MATRIX  OF  OBSERVATIONS  CXX  STORED  IN 

ARRAY  AA 

0 C»0SS:CVARIANC£  matrix  of  signal  ANO  observation  CSX 
STORED  IN  array  39 

1 AUTOCOVARIANCE  “ATRIX  OF  SIGNS  CSS  STORED  IN  ARRAY  AA 
MA*M9  DIMENSION  RAHA.METErS  OF  AA  AND  3o  AS  DEFINED  IN  CALLING  PROGRAM 

implicit  real  «8(A-H.0-Yt,L0GICAL(Ll 

COMMON  /TICa/  3?(12)tCPU0).OP(12)»ail6).CIlS)tOtlSlTF3ll6l,FC(l61 
l.TRANSil2)tSRl5l,CD,CaOfVvM,SMI,RAN.ZI5tA0OOltKRP(10I.SRSP(i2l,IOP 
2112*7)  ,K,M  lb  I fMRKl  16  I t ICI  IBf  71  ,,NN  tMM.NM.NNN,  I M ( 1 2 ) , NN  ^ . NV  * M * , I CM  . 
3NJ,MVV.MMV,I01 . ITZ.IFl 

COMMON  /C.rCCV/CI  (121  ,CRl  SI  ) ,SIG.MAO(  300  ) , S I GM  A 1 300  ) , K I C'SIfNl.LOCAL 
dimension  AAlMA,MA).3Bl.M3fMA) 

DATA  GM,R£/3. 33014. 6371. 003/ 

Input  point  p 


JV 

- NJ*“V 

IF 

(MOLD)  220,221,221 

220 

I = 

1 

CO 

TO  222 

221 

I = 

NJ»  1 

222 

ILL 

= 0 

IJ 

= 0 

214 

IJ 

- I J . I L L 

ILL 

= D 

DO  229  J=l,7 

IF  IILiPd.JI)  229,229,229 
229  ILL  = ILL.I 
228  CONTINJc 
RB  R BP( I I 
RC  = CP  I I ) 

CR( 10)  ^ GM/ ( RE»OP (1)1 **2 
RFB  = PB(  I ) 

RFC  = FC ( I ) 

CR(2)  = DP(  I ) 

J R 0 
XT  = 0 

212  IJ  = lJ*<r*ILL 
TJ  = J 
ITT  > 0 

RBR  = Ra*TJRRF9 


I 

1 


1 


non  non 
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RCR  = RC*rjsRFC 
CR(*rl  = DSIN  IR&R) 

CR(6)  = JC3s  tRjRl 
If  (MJLJI  253. ^5'., 253 

In°UT  point  3 

251  II  = 1 

IIJJ  = 0 
CO  TO  213 
253  II  = I 

IIJJ  = IJ 
210  R33  = 3P(  I II 
RCC  = CM  I I I 

Cr.  ( 1 11  = &“/(  RE.OP  ( ! I ) 1 »-2 
RF3U  = F 3(  1 1 I 
RFCv.  = FCI  II  I 
CR(3)  = 3P(II) 

IF  ( I-III  223,22t,223 
223  Jj  = J 

CO  TO  203 
22<»  Jj  = j 
206  Tjj  = JJ 

IIJJ  = IIJJ.ITT-JLL 
JLL  = 0 
DO  230  J'<=1.7 

IF  t IJP(  II. JO  I 231.230.231 
231  JLL  = JLL»1 


CG^PUTATIjN  of  SPrlPRICAL  DISTANCE 


233  CONTINJE 

RS6R  = ^3o»T J^PPPSe 
RCCR  = RCC*IJ-svFCr, 

SS  = RCCR-RCR 
CRlal  = 3SIN  (R3?R) 

CRI7I  p 3C0S  tP32R) 

CRI3)  = 3SI'.  I SS) 

CRC5I  = 3C3S  (SSI 

CP.ll)  = CR  (41  PCR  ( 5 ) *CR1  6 IPCR  ( Y ) *CR(  9) 

covariances 

IK  p 1 

304  KI  (6  I = rL'0(  I . IK  I 

IF  (KI(ol)  303.313.3:2 
300  (.RITE  (e.  333)  I.Kl  ( •- I 

303  FORPAT  (/JOX.'EiPJR  IN  INPUT  SPECIFICATION  POR  PROFILE*.  13. 

302  HI  = Ij*iK 

IF  (I-III  319.320. 319 
323  IF  ( J-JJI  319. 321, 319 
321  IKK  = IK 
GO  TO  3u7 
319  IKK  = 1 

309  X I ( 7 1 = I:P(  H . IKK  I 

IF  (<I ( 7) 1 3:5.3:7.33d 

30  5 WRITE  ( 5. j3  5 I I I .k I ( 7 I 

303  JJJ  = IIJJ».K\ 

CALL  CCV3A 

CALL  CO.CX  ICOV) 

IF  I JABSI  CO.  I 1 I 53’.oj3.5j<* 

533  COV  = 


nno  n ononor»noo  ooo 
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63<i  IF  IHOLO)  315,317f316 

316  AAdlItjJJ)  = CaV 
GO  TO  307 

317  BadlItJJJl  = CQV 
SETUP  OF  CC^^PLETE  “aTRIX 


307 

IKX 

= IX<»1 

IF 

dXR-71  309,309,310 

310 

IX 

= IX,l 

IF 

dX-T)  304,304,305 

305 

JJ 

= JJ»1 

ITT 

» 1 

IF 

tJJ-XXPdld  2C3,20 

209 

II 

= II*l 

IF 

IMQLOI  252,252,253 

253 

IF 

III-JV)  213,210,25*, 

254 

GO 

TD  211 

252 

IF 

(II-.NJI  2 10,210,211 

211 

J = 

J*  1 

XT 

= 1 

IF 

(J-XXP(Id  212,213,, 

213 

I - 

IF 

(MULO)  225,226,226 

225 

IF 

d-\JI  214,214,215 

215 

GO 

TO  227 

226 

IF 

(I-JV)  214,214,227 

227 

RETURN 

END 

SU9R0UTI*J£  ORDER  ( AA  t M A t NDEC  I 

SHIFTS  THE  U°PER  TRIANGULA®  “ART  0®  THE  REAL  SYMMETRIC  MATRIX  AA 
INTO  THE  FIRST  .MA6IMA-D/2  STORAGE  LOCATIONS  OF  THIS  ARRAY  ANLi 
VICE  versa 

**A  DIMENSION  OF  AA  AS  DEFINED  IN  CALLING  PROGRAM 

M ACTUAL  DIMENSION  OF  MATRIX 

NOEC  -I  TWO-DIMENSIONAL  ARRAY  INTO  ONE-D IME NS  I ON AL  ARRAY 

I OINEJIMENSIONAL  ARRAY  INTO  ThO-OI''£MSI  ONAL  ARRAY 

implic  it  real  s>9  I A-h,0-Y  ) tLOGICAL  (L  I 
DIMENSION  AA(MA,ma) 


IF  INOECI  300,300,301 

TWO-DIMENSIONAL  ARRAY  INTO  uNe-DI MENS lONAL  ARRAY 


300  I = 1 
J = I 
K = 1 
N = I 

304  AAd,Jl  > AA(R,NI 
X = <*l 

IF  IK-NI  302,302,303 
302  I = I»l 

IF  ( I-MA)  304,304,  305 

305  J = J»l 
I * I 

CO  TO  304 


onn  oonooooooonn  r»oo 


303  N = N*l 
K = 1 

IF  IN-MI  302,302,307 

307  CO  TO  308 

ONS-OIlcSSIJNAL  4r.57Y  I.NTQ  T',J0-0  t ME  SS  I CNAL  ARRAY  * 

301  mm  = Sf(M,l)/2 
MC  = MM/MA 
MR  = M«-MC*MA 
I = M 
J = “ 

K = MR 
N = «C*l 

IF  tR)  309,310,311 

309  WRITE  (0,312)  'J,K  i 

312  FORMAT  (lOX, 'ERROR  IN  ROSI T I ON ’ , 2 1 3 ) = 

310  R = Ma 
N = N-1 

311AA(I,w)=AA(is,iN)  * 

K.  » <-l 

IF  (K)  313,313,319 
319  I = I-l 

IF  (I)  315,315,311 
315  J = J-l 
I = J 
GO  TO  311 

313  S = N-1 
K = MA 

IF  (*n  308,308,319 

308  RETURN 
END 


subroutine  oats  (NRL,0N,G0,G0D, I Z I) 

READS  TmE  SI"ULAT£0  'cRRORLESi*  OaTA  FRO'*  A DIRECT  ACCESS  FILE  AND 
CORaJRTS  THcSE  VAi.  J£3  3Y  OIFRERcNT  RINGS  OF  ERR-RS.  IN  I HE  PREScNT 
SETUP  SIX  “EaSUaEMCnTj  aRc  aEaO  FOR  EACH  POINT,  ThE  FIRST  T-iRt  = 

OF  aHICH  are  used* 

RRL  .NJM3ER  OF  REOO.-OS  TO  3E  REaO  =RCM  EAuh  FLIGHT  FILE 

izi  -1  normal  deviates  «!Th  variances  gn,go,coo  are  aodeo  '0 

TRUE  VAL'JEj 

0 YULE  Tl'*£  SERIES  .liTrt  CORRELATIONS  0.733  A‘,0  0.307  IS 
ADDED  TO  GRAVITY  A.,0MALIES 

1 LINEAR  SYSTEMATIC  cRROR  IS  ADOEO  Tl,  GRAVITY  ANuMALIES 

implicit  real  r8(a-h,o-y) , logical (L I 

common  /t;Ca/  3P(12).CP(12),0F(1Z). 8(151, oil  51. 0(151, FR(16).FC(15) 
1,TRANS  ( 12  I .SR  ( 5 I :,0,  VVM,  Sm  t an  ,2(  5,  ADOO  I 1 1 2 I ^5P  . 12  1 , I2P 

2(12, 71, 9<(l3l,XA<l 151.10(15.71 .NN ,MM,NM, NNN, I N I ( 1 2 ! . NNV . NV . MV , I C< , 
3NJ.MVV.-Mv , 131 , I t; , 1= I 
COMMON  /ShT/  z;( 3,51 
OI—ENj  to  ( .(  lOl  ,..3[  101  .>juD(  10  1 

read  data 

DEFINt  file  11(7. 2900, L. <11 i 
0£RINE  FUE  1 2 ( T . ^ -GO  , L , < . 2 1 
OEFi'.t  file  1 Ji  r..--G-.L.<l  )1 


■k 
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502 


503 


501 

503 


DEFINE  file  ltl7.2HOO,LfK14) 
DEFINE  file  15(  7t2‘tO0,LfKl5l 
DEFINE  file  1 1, 1 7 t <;  ^00  , L 1 6 I 
DEFINE  file  I7(7»2<.00,LfH17) 
define  file  l8(7,2^00iLtM3) 
DEFINE  file  1 9 ( 7 . 2 NOO * L t < 1 9 » 
DEFINE  file  20 ( 7 . 2 ‘►00 1 L t <2 3 1 
DEFINE  file  21(7  .2‘.30.Li<2l  I 
U = 591 
IKK  = K.KL»600 
AM  = 3. 

I 3 I 


IX  = KKK( 1 1 
J1  = 1 

JK  = TUANS ( I ) ♦10. 

JJ  = 1 

J2  = Jl*599 

READ  (JK*  JJ>  ( 2 ( I «J I f J: J1 t J2) 
J1  = Jl»500 
JJ  = JJ»l 

IF  (JJ-K<L)  SOOfSOOtSOl 
I = :♦! 

IF  (I-NJ)  502,*i02.503 
IF  ( I2II  507.508.509 


c normal  deviates 
C 

507  DO  50o  i=l.NJ 
IK  =■  fa»KK(  I I 

K » -1 

JK  » TRANS! 1 ) 

DO  505  J = 1.K.6 
K = K*l 
S = G3(JK) 

CALL  GAUSS  ( U.S.AM.vi 
22(1.11  : 2( I . J > »V 
S = G00( JK  ! 

CALL  GAUSS  ( I2.S.A«.V) 

22(1.21  = 2(I.J*l)»v 
Call  gauss  ( i2.s.a9,vi 
22(1.3)  = 2(I.J*2)»v 
DO  535  JJ=1.IX 
II  = K«IX.JJ 

505  2(1.11)  = 22( I . JJ) 

506  continue 
GO  73  513 

C 

C CORRELATED  ERRORS 
C 

508  DO  511  I’l.NJ 
IK  = 6SKK( I ) 

K = -1 

JK  = TRANS  ( I ) 

S ::  GD(JK) 

CALL  GAUSS  ( I2.S.AM.U) 

CALL  GAUSS  (:2.S.AM,UU) 

V = -1.1 
VV  = .5 

DO  512  J = 1.K.6 
K = K*  I 

CALL  YULE  (U.UU. V. VV . S. SS. I 2) 
22(1.1)  » 2( I.J)»SS 


-f-%* 


o o o 


zzii,2)  = Z( r,j*i) 

ZZ( I t3 ) = :( I , J»2) 
DO  512  JJ=l,IX 
II  = K‘IX»Jj 

512  Zt I. II ) = ZZ( I, JJI 
511  CCNTI.MUc 

SYSTE>^ATIC  ERKDRS 

GO  TO  510 

509  00  513  1=1. S3 
IK  = 5-K<l 1 1 

K = -1 

JK  -=  TRaSs  I I I 
S = 5. 

SS  = K<L>100 

s = s/ss 

SS  = 0. 

DO  51-.  J = 1 . 1 < ,6 
K = K»1 

SS  = ss*s 

ZZ( I.l ) = Zl I. Jl »5S 
ZZ( I .2  I = Z(  I , I I 
ZZ( I ,3)  = :i I , J*2) 
00  5l<.  J J=  1.  U 
II  = K»IX*JJ 
51«.  Z(  I,  II  1 = Zi(  I .JJ) 

513  CO'.TINJc 

510  00  515  1 = 1. SJ 

IK  = KK  1 I ) 4=KKK(  I ) *l 
00  5’t>  j=:k.ikk 
516  Z( I. Jl  = 9R939. 

515  CO'iri'.Je 
R £ "UAN 
ENO 


SoaROJTriE  SCA9  IAA.9A.  I''li.V'*.RR.  ICKEl 


COMPUTES  0'l£  profile  C=  POINT  OR  MEAN  VALUES 
Aa  PREOICTION  -'ATRIX 

1.''“.  sequential  N.J“jEa  0=  VLAN  OR  POINT  VALuE 

RR  -RRAr  TO  STOr^E  TmE  C0'*PUTEJ  vacUES  a NO  T‘“?IP  COORI 

ALL  OTREP  OUAnTIT.Es  A.Pl  AS  JE^I'IEO  I N .IAlLINO  PPuORA- 


'.nates 


t«PLt::T  ACM.  -3(  A— I,  ) .logic  AL  (L  ) 

CO-NON  /TIC-/  3P  ( 1 0<  .C  = I 12  1 .OP  ( I ’ I . 3U6I  .L  ( 15  I .01  1<J  1 .F->  I Is  I 
1 . TPaNS  ( I 2 1 , j 3 ( 5 I ...  0.  COO.  VV*.  S'*  I . P AN  . ; 1 5.  a : 3..  I .'.N?  < 1 2 I . V.\<P  I 
Z(  12.71  .■■<  1 l->  1 .»  K ' 1 I 5 I . I Ji  15.7)  ...N  .RM.NR.  NNN.  INI  I I 2 I .N.NV  . Sv  . 
3NJ.  “Vv  . . I - I . I . I F I 

OrVNSION  AhMA.aa  I.RRIoIZ.o) 


. = C I I 
1 0 I . I ’’  P 
■■  V . I c < . 


CO'^PASE  initial  “OInTj  OF  FLIGHT  ANJ  I NT  ER  POL  A T I ON  pPOFILcS 


JJKK  = I Xr  I 1 I - 1 I =<  KN  I 1 I - 0« '«HV»''VV 
IF  lion  S22.s2i.i23 
623  SV**  = V- 
GO  TO  52 A 
622  SVP  = 0. 

62  5 ScH  = •'V M'*V/2  I 
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BLO  = C( rCKfcj-Svx 
BLAT  = BlIC^E) 

BL  = BLO-SU'i*  VV''/I  OC33  ( BLaT  ) )-S»'</2. 

KH  = KKt IC<t I 
00  616  1=1 tNJ 
616  INI  I I ) = 0 
IM  = 1 

607  SUM  = IM-1 
SUMI  = 0. 

BLO  = 3L0*SVM 
KR  = NNV 
SU‘’S  = aR 
IR  = 1 
3L  = 3L*SvM 
NC  = TRANSm 
CPR  = Cl  NCI 

CCD  = t BL-CPO  I <=DCJSt  3LAT  !/VVM 
IF  ICJD)  600tCs01t602 

600  IF  (C00».5l  603,601,601 

603  NC  = COD-. 5 
KRR  = iCK»NC 

IF  (KHRJ  6JA, 604,605 

604  IC  = ICne-NN 

WRITE  16.O06I  IC,IM,Imm 

606  FORMAT  ( // 2a  , • I N T6  RF  OL  AT  I ON  PROFILE  N0',:4.'  °QI.NT  NO',!**,' 

ISEOUcNTiAL  NT,I4,/2X,  'MAS  ■(OT  BEEN  CJ‘'POTE0  BECAUSE  NO  '^EASURE.M 
2ENTS  ARE  AVAILABLE’,/) 

RRil'-M,!)  = BlO 
RRIIMM.Z)  = 999999999. 

RRIIKM,3I  = 0. 

IM  = IM.l 
IMM  = IMM*1 

IF  IIM-RM)  607,607,630 
630  GO  10  615 
6C5  SU“  = -i.C 
COO  = XK 

CN  = lOG.-SOMFlOO./COO 
RRl IMM, 3)  = Cn 
SU^S  = KR»NC 
IR  = IR-NC 
IC  = ICRE-NN 

WRITE  (6.603)  IC  , I “,  IMM.C.N 

608  FORMAT  I // 2X  , • I MTE  RP  OL  AT  I ON  PROFILE  NO  ' , I 4 , ' POINT  NO’,:-,,' 
ISEOUE.NTIAL  N0’,I4,/2X,  ’mas  aEEN  CuMPJTED  aITh  ONLr’.FE.l,'  P£RC 
2ENT  OF  TWE  -1 E 4 SURE '’L  NT  S RECUIREO',/) 

GO  TO  601 

602  IF  1C0D-.5)  601,601,609 

609  DO  olO  1=1, NJ 

COO  = I &L-CPP ) «OCOSt 3P( I ) )/VVM 
NC  = LOO*. 5 

610  IM  ( I ) = NCPRKXP  ( I ) 

601  IF  t lOI ) 625,025,626 
C 

C MEAN  values 
C 

626  CALL  FIX  I 4a,MA,SoM, JJ) 

Sji:  = SJMI.SU" 

IF  1 INI  I 1 )-JJRR)  ol7,ol3,6l3 
613  SUMS  = IR 
COO  = VR 

IF  I IR-aR ) 620.621 ,620 
620  CN  = SUMS* 100. /COO 


O o o o 
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621 

617 

612 

619 


C 

C 

C 


629 


625 


62  7 

628 

631 

632 
615 


RK( = CN 
IC  = 

WRITE  16,533)  1C,I“, 

I.M  = 

30  TO  619 

IR  = IR»1 

DO  612  I=1,NJ 

INK  I ) = I 91  (1  I ,K^<R  ( I ) 

IF  I IR-sR)  626,625,019 

SU'^I  = SJII/SJ'-S 

RRI  I««1,;  ) = 2l3 

RR(  I >••■),  2 I = S J-I 

IM'1  = IM''»1 

I.M  = I'*.  I 

IF  ( I.'1-<M)  677,607,629 
30  TO  615 

POI.ST  VALOEi 

Call  Fix  1 AA ,«A, S'j-,  j j 1 

Rk(IMm,ii  = 

RR(I“m,,:)  = S^M 

3l3  = SlO-VW/oCOSOlAT) 

Imm  = I“''»l 

IM  = I“«l 

DO  6.’ 7 t = l,9J 

IM  I I ) = I 91  I 1 ) 1 I) 

IF  ^IM-^'')  623,628  ,01  5 

IF  1 19  I t 1 ) ♦.s.I.N  P ( 1 I - J JKR)  6 2 5, 631 , 63  1 

00  632  : = 

RRIl,2)  = 999999999. 

return 

EN3 


9 


SL)?R0oTri6  FIX  I AA  , “A  , 3 J**,  J J ) 

pE9FM9“S  one  CPCRAriJ'.  QF  T«F  PREOrCTlCN  matrix  aa  09  ’’.-lE  '265cPVAnw9 
field  2.  n.£5JL”  IS  jTOREO  1.9  jOm, 

i.mplicit  real  *0  (a-o ,0-x ), logical  (L  I 

COM.MO.9  / T ; Cx  / 2P  ( 1 2 , 2 P 1 I : I , DM  ( I 2 ) . 3 ( 1 . C I 1 6 ) , C 1 I I , P 5 ( 1 1 ) , c ; [ I ‘ I 

l,TRA9S(l2),'SRl5J,,.2,CLi.7.V'.M,3Mi,6A-|,Z|5.-2'j3).<'‘P(12),-.  \ MP(!2i,l 
2(  12,7;  ,X9t  UJ  I lo  I , rC(  lo,  7)  , 99  , 9,.9  , I . I ( I 2 I , 9 9V  , 9'.  , -•V  , ' . X , 

39J,MVV ,“Mv. I 31  , I T1  , ; F I 
3I''t93I09  aa(ma,“a) 

1 = 1 

733  SUM  5 3, 

J = 0 
II  = 1 
JLL  = 3 
701  ILL  = 7 

KRX.X  = .X<.\P)II) 

XO  = M9\,pxkk< 

ILL  = ILl.-XO' I9l  I I I ) 

JLL  = JLL-X.'.XX 
30  7 00  !L  = 1,M', , 

ILL  = lLL*Xo 
JLL  = JLL*<\Xx 
30  700  i:  = l,X,xK\ 
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I 


j = JLL*rZ 
JJ  - ILL»IZ 

700  SUM  = SuM»4AI  IfJ)<=Z(  Ilf  JJ) 

JLL  = JH.*KM<H 

II  = Ilfl  ■ - ' 

IF  (II-\J)  791r701,702  ■< 

702  SRI  I)  = SUM 
I » 1*1 

IF  I r-My)  703,703,  70<> 

70‘>  SUM  I 0. 

00  70u  1=1, MV 

70S  SUM  r SUMfSfU  I)  ■ 

SUMI  = r/ 

SU“  = SuVSu“i 

RfcTUR\  ' 

ENO  i 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 


SuaROUTINE  C0.MPA  I MR  , RE  , SUM  , SUMS  ) 

CCTERMINES  THE  OIF-FERENCES  eETHEEN  TRUE  A.V?  ESTIMATED  VALUES  OF  THE 
SIMULATES  mooel.  Tn-  SUM  OF  SJUaRES  AHO  The  STA.'vOaRO  OSvIATIOV  ARE 
STORES  IN  SUM  A.NO  ShMS 

RR  array  containing  the  TRUE  ANO  THE  ESTIMATED  VALUES 

RE  total  NUM3cR  OF  ESTIMATEO  VALUES 

ITZ  total  numser  of  True  values 

SuM  SUM  of  scuares  of  differences 

Sums  variance  of  differences 

IMPLICIT  REAL  I A -h  , 0- Y ) , LUG  I C AL  I L ) 

COMMON  /TtCR/  ?P (IZ!.C0|12I, OP (121, 3(161, CilS), oils), FOIlbl.FCIlEl 
IfTRANSI 12) ,SRI 51 ,CO,COO, Vv",SHI,RaN,Z( 5,aZ0D) ,RrP( 12) ,Y\<PI 12i , nf> 
2ll2,7),R<(lo),KRRIl6),tOIIo,7) ,nN ,mm,NM, NNN, I N I 11 2 I , NNV , NV , m V , I C « , 
3NJ,MVV,mMv , IOI . i TI , if  I 

dimension  RRtol2,t) 

SUM  = 0. 

RKKR  = 0 
IDO  = 3 
ILL  = I 
J = I 
JJ  = 1 
NNNN  = NNN 
706  ELAT  = 3INN.NN) 

702  IF  I CA&SI jLA T-RR I J,A ) 1-.31 I 701,701,700 

700  J = J*1 

IF  IJ-IFIl  732,702,733 

703  IF  liLL)  71„,71h,713 
713  write  Ii,710) 

710  FORMAT  I ///2 < , 'TSUF  AND  ESTIMATED  VALUES  CO  NOT  HAVE  COMPATIRlE  CO 
lOkDINATtS* ,/// ) 

GO  TO  709 

701  BL3  = RRIJJ.n 

711  IF  ( OaPS(.(Lj-RR(  J,  5)  ) -.01 ) 709,70m,  712 
709  KK<N  = RSINNNN! 

KKRK  = K'RKtRMSN 
DO  70S  I=jJ,K.Y.M9 
II  = J-JJ»I 
RRI I ,6  ) = RR I I I ,ol 
70S  RRI I ,R ) = RRI  I I , 5) 

ILL  = -I 


ooor>oooooor'»o  oooooooo 


li 


\ 

i 

1 


51 


JJ  = 

J : J*KK<.'.I 

NNNN  = NN.N'Ul 

IF  (NNNN-N'I)  70tj,70b.707 

707  IF  (ILL)  7 !<..  7 I'f  ,7  13 
714  00  703  i = l ,i<e 

RK( I f 4 ) = RR ( ! I 3 I 

RR(  I .5  ) = 3^  ( I t6  I 

IF  IRRII.J)  -999=)'.i97  99.  ) 716,717,716 

716  SS  = RRl  1 , 21  -N=,  I 1 , 51 
SOM  = SU’1»SS5S3 
RRM,6)  = S3 

GO  rO  703 

717  IDO  = lOO-l 

rr(i,6)  = 999999997. 

RM  1 ,3  1 = 0. 

708  CGNURJc 

S3  = RS-lOn 
RE  = SS 

SU-S  = OSORT  (SUM/SSl 
GU  TO  709 


712 

J A 

J*  1 

IF 

1 J-  IF 

I)  711,711.715 

715 

IF 

1 I L L 1 

709,709,713 

709 

RETURN 

END 

SUSROjTINc  TLILE  (U,U'J,V,VV,S,S3,IZ1 

CORRUTES  - STaTI0\4RY  TI“E  SERIES  '.ISING  YULE'S  SCmE-E 

u.'ju  coefficients  of  The  series  for  The  T«0  pRECEDINO  ECEnTS 

V,VV  the  two  RRECEOiNG  EVENTS  OF  I HE  SERIES 

S VARIANCE  OF  the  OAUSSIAN  PROCESS  U3£o  TO  GENERATE  THE  SEhIES 

U RA,DC'*  NuHEE^  TO  STA^T  GAUSSIAN  PROCESS 

SS  transfer  of  C'O.HPUTEO  E/E.NT 


AM  = 3. 

CALL  OajSS  ( i:,3.AM,6) 
SS  = E-uPV'uUPVV 

UU  = w 

0 p SS 
RETURN 
END 


SU3R00TINE  CO.AA 

A listing  of  This  SJ.»RCjTINE  can  re  F"'u*»u  in  TSC^p-vN1*JG  11/761 
OSU  RlPl’RT  no  206.  The  C'’k  = ECrt^NS  .'-'.JCO  RE  HAOE 

(PERSONAL  C-HhuM:  AT  tuN  3Y  C . C . I 3 CH  ER  N I N 0 i : 

page  59. line  23  i --R  L 1 I -•  R ^ • I 0 . > : 5 . • ( 6 . ♦ 00  *R  L I I «RL  1 I * R L I I RR  N I aR  M -'R  N I A F , L 
PAGc  sO, statement  r;  Ch-'-,Gc  3I.N 
PaGF  6 3 . S T a T E " r N T 13  ChA.Gc  S I 'SN 


THE  page  .UMRERS  refer  TO  AqOVE  REPoOT 


Appendix  B:  Sample  Computations 


K M 

p L e 

INPUT 

99^00000*00 

.9710^003* 

7 ■ 

‘.9b. 

6 10 

7 • 

60. 

10.  0. 

a. 

40. 

10.  0, 

9. 

40. 

10.  0. 

0. 

60. 

10.  0. 

I • 

60. 

10.  0. 

2. 

60. 

10.  0. 

3. 

60. 

10.  0. 

1 

2 10 

150 

7.5 

t 

o 

0.  15 

8.5 

60.  5 

0.  10 

PAHAMElfHS  fPrCIEYINi'.  IHK  hf  CUE  f - VA  R 1 ANC  F MCOEL 


INVfoSITN  OF  C0V4RIANCF  HAlRIX  OKAY 


Sllf*  OF  I0U4PF  S:  ARA  .«  3 

SIANDARD  nrvlATION  FPt.O  ?5  M <•  * c ,j(,  f f,p  m n . 4,^,; 

FORMAL  SIANOARO  DFVIA1ION  FROM  C OL  LOCA  T I ON  : 6. OS  Mf.AL 


